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Abstract 

In  this  project,  methods  of  model  reduction  that  integrate  feedback  active  flow  control  with 
applications  to  nonlinear  convection  and  turbulent  flows  governed  by  Navier-Stokes  equations 
are  developed.  A  new  methodology  which  extracts  boundary  conditions  in  reduced  order  proper 
orthogonal  decomposition  (POD)  and  finite  difference  models  is  developed.  A  new  model 
reduction  method  based  on  empirical  data  and  balanced  truncation  was  developed  and  applied  to 
nonlinear  Galerkin  models.  Based  on  this  method  a  new  empirical  Hankel  norm  model  reduction 
algorithm  is  proposed.  These  methods  are  applied  to  a  prototype  nonlinear  convective  problem 
governed  by  the  two-dimensional  (2D)  Burgers’  equation.  The  reduced  models  are  used  in  the 
design  of  robust  boundary  controllers  that  achieve  tracking,  and  implemented  on  the  full  order 
Computational  Fluid  Dynamics  (CFD)  models. 

POD  and  balanced  truncation  are  shown  to  be  optimal  in  the  sense  of  distance  minimizations  in 
spaces  of  Hilbert-Schmidt  (HS)  operators.  POD  has  been  shown  to  be  optimal  in  a  broader  sense 
than  is  reported  in  the  literature.  The  optimality  is  in  the  sense  of  distance  minimization  in  a 
space  of  integral  operators  under  the  Hilbert-Schmidt  norm.  A  connection  with  balanced 
truncation  is  found.  In  particular,  balanced  truncation  is  shown  to  be  optimal  in  the  sense  of 
distance  minimization  albeit  in  a  different  space  of  integral  operators  under  the  Hilbert-Schmidt 
nonn  when  the  so-called  ’Curtain-Glover’  balanced  realization  is  used.  This  is  a  novel  discovery 
as  balanced  truncation  is  usually  known  to  be  not  optimal  in  any  sense. 

An  algorithm  that  combines  the  extended  Kalman  filter  and  expectation  maximization  to 
estimate  the  model  coefficients  from  particle  image  velocimetry  (PIV)  measurements  for 
turbulent  2D  flows  is  developed.  The  algorithm  is  recursive  and  convenient  for  on-line 
implementation.  The  method  is  applied  to  a  2D  flow  control  problem  over  the  NACA  4412 
airfoil  using  experimental  data  obtained  from  Prof.  Glauser’s  flow  control  group  at  the 
University  of  Syracuse.  The  motivation  for  this  problem  is  minimization  of  aero-optic  distortion 
in  a  bluff-body  flow. 

POD  and  balanced  truncation  are  shown  to  minimize  different  n-widths  of  partial  differential 
equation  solutions  including  the  Kolmogorov,  Gelfand,  linear  and  Bernstein  n-widths.  The  n- 
widths  are  notions  from  metric  complexity  theory.  They  quantify  inherent  and  representation 
errors  due  to  lack  of  data  and  loss  of  information. 
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Most  model  reduction  algorithms  assume  linear  models  and  fail  when  applied  to  nonlinear  high 
dimensional  systems,  in  particular,  fluid  flow  problems  with  elevated  Reynolds  numbers.  For 
example,  POD  fails  to  capture  the  nonlinear  degrees  of  freedom  in  these  systems,  since  it 
assumes  that  data  belong  to  a  linear  space  and  therefore  relies  on  the  Euclidean  distance  as  the 
metric  to  minimize.  However,  snapshots  generated  by  nonlinear  PDEs  belong  to  manifolds  for 
which  the  geodesics,  when  they  exist,  do  not  correspond  in  general  to  the  Euclidean  distance.  A 
geodesic  is  a  curve  that  is  locally  the  shortest  path  between  two  points.  We  propose  a  model 
reduction  method  which  generalizes  POD  to  nonlinear  fluid  flows  corresponding  to  manifolds 
which  have  a  differentiable  structure  at  each  of  their  points.  The  algorithm  is  applied  to  a 
prototype  nonlinear  turbulent  flow  on  the  Wortmann  FX  60-100  airfoil  governed  by  the  Navier- 
Stokes  equation. 

I.  Introduction 

Recent  advances  in  actuators,  sensors,  simulation,  and  experimental  diagnostics  bring 
applications  such  as  suppression  of  acoustic  tones  in  cavities,  separation  control  for  high  lift,  and 
trajectory  control  without  moving  hinged  surfaces  within  reach.  However,  these  applications 
require  the  integration  of  feedback  control  because  of  the  need  for  robustness  to  flight  condition 
and  vehicle  attitude,  precision  tracking,  overcoming  low-fidelity  models,  or  moving  a  system 
away  from  a  stable  solution  or  limit  cycle  as  efficiently  as  possible.  Feedback  control  strategies 
in  which  the  bandwidth  of  the  controller  is  commensurate  with  the  time  scales  of  fluid  flows 
are  attractive  because  they  offer  the  possibility  of  improved  perfonnance  and  reduced  control 
power  required  through  control  of  unstable  structures  in  the  flow  field.  Unfortunately,  models 
that  capture  the  relevant  dynamics  of  the  input-output  system  and  are  amenable  to  control  design 
are  difficult  to  develop  [1]. 

There  has  been  significant  interest  in  model  reduction  for  the  purpose  of  control  design 
[2][3][4][5][6][7][8][9][10].  One  such  application  of  reduced  order  modeling  is  control  design  in 
the  context  of  aerodynamic  flow.  Aerodynamic  flow  control  is  a  research  area  of  great  interest  to 
the  Air  Force  and  the  fluid  mechanics  community.  Presently  considerable  research  efforts  are 
working  with  feedback  control  law  design  for  systems  described  by  PDEs  that  need  a  very  large 
number  of  states  to  accurately  simulate  their  characteristics.  However,  recent  advances  in  the 
design  of  actuators  and  sensors  can  be  leveraged  for  better  system  control  only  if  the  control 
design  methods  provide  a  reliable  low  order  controller  [11].  Additionally,  simulation,  and 
experimental  diagnostics  are  making  applications  such  as  the  suppression  of  acoustic  tones  in 
cavities,  separation  control  for  high  lift,  and  trajectory  control  without  the  need  to  move  hinged 
surfaces  a  possibility  [12].  However,  these  applications  require  the  integration  of  feedback 
control  because  of  the  need  for  robustness  to  flight  condition  and  vehicle  attitude,  precision 
tracking  overcoming  low-fidelity  models,  or  moving  a  system  away  from  a  stable  solution  or 
limit  cycle  as  efficiently  as  possible  [12]. 

Unfortunately,  it  is  very  difficult  to  create  models  that  capture  the  relevant  dynamics  of  the 
input-output  system.  For  example,  computational  fluid  dynamics  simulations  can  provide  good 
solutions  to  a  discretized  version  of  the  Navier-Stokes  equation  [2],  [13].  However,  accurate 
simulations  for  simple  shapes  such  as  two-dimensional  airfoils,  or  complex  shapes,  such  as  a  full 
vehicle,  require  several  thousands  to  millions  of  states.  Therefore,  the  simulation  results  are  not 
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directly  useful  for  control  design  [12].  Complexity  in  the  model  is  a  legitimate  need.  The  large 
number  of  states  is  necessary  to  capture  important  flow  features  that  occur  at  extremely  small 
spatial  scales.  Although  these  small  flow  features  might  seem  insignificant,  if  they  are  not 
captured,  it  is  not  possible  to  analyze  if  they  are  necessary  in  securing  the  closed  loop  system’s 
overall  stability  [11]. 

POD  has  been  extensively  investigated  in  distributed  parameters  systems  due  to  its  order 
reduction  capability  [8]  [9][  1 0]  [  14]  [  1 5],  and  balanced  truncation,  which  is  a  simple  yet  efficient 
model  reduction  technique  widely  used  in  reducing  model  orders  of  high  order  linear  systems 
[16]  [18]  [19].  POD  models  of  only  a  few  dozen  states  can  often  accurately  capture  the  input- 
output  behavior  of  systems  that  have  full  order  system  models  of  thousands  of  states  [12].  In 
addition  to  using  the  POD  method  in  conjunction  with  model  reduction  techniques,  the  idea  of 
using  empirical  gramians  is  growing  in  popularity  for  use  in  an  approximate  balanced  truncation 
[20] [5]  [21] [17].  Further,  some  work  has  been  done  on  finding  nonlinear  empirical  gramians  for 
balanced  truncation  [22] [23]. 

In  fluid  flow  configurations  it  is  not  uncommon  for  discretized  flow  models  to  describe 
thousands  to  millions  of  state  variables,  for  example  if  one  uses  a  linear  quadratic  regulator 
(LQR)  control  formulation,  roughly  10  "  Riccati  unknowns  need  to  be  calculated  for  a  discretized 
flow  model  describing  106  states.  Existing  computing  power  and  computational  algorithms  are 
not  capable  of  solving  an  LQR  problem  of  such  large  dimension.  For  dynamical  models  that  are 
very  large  scale,  such  as  those  describing  fluid  flow  configurations,  it  is  apparent  that  the  order 
of  the  system  must  be  reduced  prior  to  control  law  design  [24].  This  prevents  us  from  using 
closed  loop  model  reduction  strategy  wherein  the  system  is  part  of  a  closed  loop  system  with  a 
controller  [25], 

In  the  area  of  fluid  mechanics  controls  must  often  be  fixed  to  the  boundary  of  the  problem 
geometry.  For  example,  control  of  flow  separation  over  an  airfoil  requires  that  actuation  and 
sensing  be  done  on  the  airfoil  surface  [12].  The  problem  geometry  used  for  this  project  is  one 
example  of  a  case  where  control  is  restricted  to  the  boundaries  by  physical  necessity. 

The  rest  of  the  report  is  organized  as  follows.  In  section  II,  POD  is  shown  to  be  optimal  in  a 
more  general  sense  to  what  is  reported  in  the  literature.  In  particular,  POD  is  shown  to  be  optimal 
in  the  sense  of  distance  minimizations  in  space  of  Hilbert-  Schmidt  or  trace-class  2  integral 
operators.  In  section  III,  balanced  truncation  is  shown  to  be  optimal  in  terms  of  approximation  by 
finite  rank  operators  in  the  Hilbert-Schmidt  nonn  when  the  Curtain-lover  balanced  realization  is 
used.  In  section  IV,  a  prototype  nonlinear  convection  problem  is  introduced.  A  new  empirical 
balanced  truncation  is  developed  in  sections  V  and  VI.  The  method  is  applied  to  nonlinear 
convection  in  section  VII.  A  new  empirical  Hankel  nonn  model  reduction  is  developed  and 
applied  to  the  nonlinear  convection  prototype  problem  in  sections  VIII  and  IX,  respectively.  In 
section  X  a  HJ  flow  controller  that  achieve  tracking  for  the  nonlinear  convective  flow  is 
designed  based  on  the  reduce  order  models  and  applied  to  the  full  order  system.  In  section  XI, 
model  estimation  and  identification  algorithms  that  combine  the  Kalman  filter  and  expectation 
maximization  for  POD  modes  of  a  turbulent  flow  over  the  NACA  4412  airfoil  are  presented.  In 
section  XII,  POD  and  balanced  truncation  are  related  to  metric  complexity  theory  with  the  notion 
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of  n-with  approximation.  In  section  XIII,  a  generalization  of  POD  to  nonlinear  PDEs  using  the 
notion  of  auto-associative  models  is  proposed.  Section  XIV  contains  concluding  remarks. 

II.  Optimality  of  the  Proper  Orthogonal  Decomposition  (POD) 

POD  has  been  used  extensively  to  determine  efficient  bases  for  dynamical  systems,  random 
processes  and  large  data  set  in  general.  It  was  introduced  in  the  context  of  turbulence  by  Lumley 
[26].  It  is  also  known  as  the  Karhunen-Loeve  decomposition,  principal  component  analysis, 
singular  systems  analysis,  and  singular  value  decomposition  [27] [28].  The  fundamental  idea 

behind  POD  is  as  follows:  Given  a  set  of  simulation/experimental  data  or  snapshots  {S^  of  a 

function  w(t,x)  ,  in  the  standard  Hilbert  space  Zr([0,r)xQ)  ,  where  xgQ  for  some  set 

(spatial  domain)  Q  of  M  p  and  T  represents  a  finite  or  infinite  time  interval.  The  n- th  POD 

vector  (4(x)  is  chosen  recursively  so  as  to  minimize  the  cost  function  [26]  [29] 

J((t>n)'-=  J  J I  $i (*, x) -  (x)  |2  dxdt  (1) 

[o,r>  n  j= i 

subject  to  the  constraints 

aJ(t)  =  ^Si(t,x)<t>J(x)dx  (2) 

Q 

J ' <f>i(x)<t>](x)dx  =  Sij,  for  i,j  =  1,2,— ,n  (3) 

Q. 

The  optimal  POD  basis  is  given  by  the  eigenfunctions  {$(x)|  of  the  averaged  autocorrelation 
function,  denoted  R(x,  x')  ,  of  the  snapshots,  that  is,  [26][27] 

f?(x,x')  =  J  Si(t,x)Si(t,x')dt  (4) 

[0,7) 

which  solves  the  eigenvalue  problem 

[  [  Sj(t,x)  Si(t,x')<f>(x')dtdx  =  A</>(x)  (5) 

n  [0.7) 

The  Hilbert  space  L~  ([0,7")xf2)  is  endowed  with  the  norm 

i 

(  \2 

\\Mt,*)\\2  ■=  Hi  w(t,x)|2  dxdt  (6) 

v[°,0  n  j 

Define  the  shortest  distance  minimization  in  the  ||  •  ||^  -norm  from  the  function  w(t,  x)  to  the 
subspace  S,  by 

H  :=inf  ||w(t,x)-x(t,x)||  (7) 

seS 

where  the  subspace  S  is  defined  as 

5  :=  jlX-  (0  (h  (x)’  ak (0  e  L2 ([°» T)) »  A  (x)  e I2  (Q),  n  integerj  (8) 

Note  that  this  distance  problem  is  posed  in  an  infinite  dimensional  space.  For  finite  dimensional 
spaces,  in  particular  for  distances  to  lower  rank  matrices  see  [30],  where  SVD  techniques  are 
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used.  To  compute  the  distance  we  view  w(t,x)as  a  Hilbert-  Schmidt  kernel  for  an  integral 
operator  T  mapping  Z2(Q)  into  L2([0,  T))  both  endowed  with  the  standard  ||  •  |, -norm,  and 
defined  by 

(t)  :=  f  w(t,  x)t/>(x)dx  (9) 

Q 

It  is  known  that  such  an  operator  is  compact  [31],  that  is,  an  operator  which  maps  bounded  sets 
into  pre-compact  sets.  The  operator  T  is  said  to  be  a  Hilbert-Schmidt  or  a  trace-class  2  operator 
[32].  Let  us  denote  the  class  of  Hilbert-Schmidt  operators  acting  from  Z2([0,r))  into  Z2( Q),  by 
C2 ,  and  the  Hilbert-Schmidt  norm  ||  •  ||  . 

Define  the  adjoint  of  Tby  T*  as  the  operator  acting  from  Z.2([(),  T  j)  into  L:(£l)  by 

<Tf  ,g  >2  '■=  f  \w(t,x)f(x)dxg(t)dt 
[o,n  n 

=  f/(x)  f  w(t,x)g(t)dtdx 

n  [0,r) 

=:  <  f,  T*  g  >, 


(10) 

(ID 


showing  that  (T*g)(t)  =  J  w(t,x)g(t)dt . 

[0  ,T) 

*  y 

Using  the  polar  representation  of  compact  operators  [32],  T  =  U(T  Ty 2 ,  where  U  is  a  partial 

*  y 

isometry  and  (T  Ty2  the  square  root  of  T,  which  is  also  a  Hilbert-Schmidt  operator,  and  admits 
a  spectral  factorization  of  the  form  [32] 

(r*r)^  =  ^A,v,®v,  (12) 

i 

where  X{  >0,  Xt  \  0  as  i  Too,  are  the  eigenvalues  of  (TT)/2 ,  and  vi  fonn  the  corresponding 

*  y 

orthonormal  sequence  of  eigenvectors,  i.e.,  (T  T)/2vj  =  X:vi,  i =1,2,  •••.Putting  Uvi  =:  cpt,  we  can 
write 

T  =  YJ^i®cpi  (13) 

i 

Both  |v,  }and  (<£>,  }  are  orthonormal  sequences  in  L2([0,  T))  and  L2(Q) ,  respectively.  The  sum 
(13)  has  either  a  finite  or  countably  infinite  number  of  terms.  The  above  representation  is  unique. 

Noting  that  the  polar  decomposition  of  T*  =  U*(TT*Y 2 ,  a  similar  argument  yields 

(7T-)Ux^®« 

i 

^  =2^, 

i 


5 


(14) 

(15) 


*  y 

which  shows  that  ai  form  an  orthononnal  sequence  of  eigenvectors  of  (TT  y2  corresponding  to 
the  eigenvalues  /L .  From  (12)  and  (15)  it  follows  that 

Tq>,=U(fT)\,=^,  (16) 

7V,=C/(7T’)5v,=V,  (17) 

We  say  that  q>{  and  v;  constitute  a  Schmidt  pair  [31].  In  terms  of  integral  operators  expressions, 
identities  (16)  and  (17)  can  be  written,  respectively,  as 

y(0  =  j"  (18) 

Q 

T 

(pi(t)  =  \w(t,x)vi(t)dt  (19) 

0 

In  terms  of  the  eigenvalues  A; 's  of  T,  its  Hilbert-Schmidt  nonn  || .  ||//v  is  given  by  [32] 

l  I 

II  rL  =fz^2]2  =  [)\\w(t,x)\  dxdt  2  (20) 

V  i  )  V  o  n  j 

Note  that  since  the  operator  T  is  Hilbert-Schmidt  the  sum  in  (20)  is  finite.  The  Hilbert-Schmidt 
nonn  is  also  induced  by  the  operator  inner  product  defined  by  (23). 

By  interpreting  each  elements  of  the  subspace  S  defined  in  (8)  as  a  Hilbert-Schmidt  operator  as 
we  did  for  w(t,x)  we  see  that  S  is  the  subspace  of  Hilbert-Schmidt  operators  of  rank  n,  i.e., 

^  =  {S  =  Y'jMWzMV-fji 0  e L\T ),  Zj(x)6f2(Q),^M}  (21) 

In  addition,  the  distance  minimization  (7)  is  then  the  minimal  distance  from  T  to  Hilbert-Schmidt 
operators  of  rank  n.  In  other  terms,  we  have 

//  =  min||r-45  (22) 

The  space  of  Hilbert-Schmidt  operators  is  in  fact  a  Hilbert  space  with  the  inner  product  [32], 
denoted  (■ ,  •),  if  A  and  B  are  two  Hilbert-Schmidt  operators  defined  on  Z2  (Q) , 

(A,  B)  :=  tr(B* A)  (23) 

where  tr  denotes  the  trace,  which  in  this  case  is  given  by  the  sum  of  the  eigenvalues  of  the 
operator  B  A  which  is  necessarily  finite  [32].  Note  that  the  inner  product  (23)  induces  the 

Hilbert-Schmidt  norm  ||^4||//s,  =  (tr^A* A^}  2 . 

In  the  case  where  A  and  B  are  integral  operators  with  kernels  A(t,x)  and  B(t,x) ,  respectively, 
the  inner  product  can  be  realized  concretely  by 

T 

(A,  B)  =  J  |  A(t,  x)B(t,  x)dxdt  (24) 

o  n 
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The  solution  to  the  distance  minimization  (22)  is  simply  given  by  the  orthogonal  projection  of  T 
onto  S.  To  compute  the  latter,  note  that  the  eigenvectors  of  ITT  j  2  and  I T*T\ 2  form 

orthonormal  bases  (by  completing  them  if  necessary)  for  Z2([0,r))  and  L1  ( Q) ,  respectively.  In 
terms  of  the  eigenvectors  v.  and  (pj  the  subspace  S  can  be  written  as 

S  =  Span  {Vj  ®  (pj ,  j  =  1, 2,  •  •  • ,  n}  (25) 

Since  the  shortest  distance  minimization  (22)  is  posed  in  a  Hilbert  space,  by  the  principle  of 
orthogonality  it  is  solved  by  the  orthogonal  projection  Ps  acting  from  C,  onto  S.  The  latter  can 

be  computed  by  first  determining  the  orthogonal  projection  P  onto  Span  { v; ,  j  =  1, 2,  •  •  • ,  n) ,  and 
the  orthogonal  projection  P(p  onto  Span  {(pp  j  =  1,2,  •••,«}.  These  projections  have  finite  rank 
and  since  the  v,  and  q>.  's  are  orthogonal  vectors  in  L2  ([0,  T ))  and  L1  (Q) ,  respectively,  it  can  be 
easily  verified  that  Pv  and  P/p  are  given  by 


(Pf)  (0  =  S  (/’  Vi  ),  Vj  (0  =  X  j /(0Vy  (0  dt 

j= 1  y=l  Vo  y 


v,.(0 


(26) 


(P,J) 


(pj) (x) = X (G’ (Pj )2  <Pj = X  J (x)Jx 

y'=i  t=i  V  n 

The  overall  orthogonal  projection  can  be  computed  as 

Ps=Pv®J 

That  is,  if  W  e  C2  has  spectral  decomposition  X  ®  v;. ,  where  e  L2(T),  ,9;  e  L2(Q)  then 

1=1 


(27) 


X;A 


X  J“l(0v/0* 

v  t=1  VO  J 


V.  ®^(^(x)^.(x)Jx)^. 


7y 


(28) 


=  J]djvj0<Pj,  3  scalers  0 

j= i 

where  the  last  finite  sum  is  obtained  thanks  to  orthogonality,  i.e.,  only  the  uj 's  and  9.  ’s  that  live 
in  the  span  of  v.  ’s  and  cpP s,  respectively,  are  retained.  For  the  orthogonality  property  we  only 
need  verify  that 

x<8)y-(Pv  ®i3p)(x®_y)  ±w®v,  xeL2([0,T)),  yeL2(Q),  u®vgS 
Computing  the  inner  product,  we  get 

<x-Px,u  >i<y-P(p,v>2=  0 

Because  P  is  the  orthogonal  projection  of  Z.2([0,  T))  onto  Span{v-,  j  =1,2,  •••,«},  and  P/f)  the 
orthogonal  projection  of  L2(Q)  onto  Span{(pj ,  j  =1,2, •••,«}  .  The  minimizing  operator  ,s0  e  S  in 
(22)  is  then  given  by 

v=-Psr  =  i>,®«  (29) 
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and 


N  1/ 
\/2 


//  =  mm|w(?,x)-^,x)|2  =  Yj  Ai  I  (3°) 

\i=n+ 1  J 

And  as  nToo,  ||r-/^7’||  \  0.  Therefore,  the  minimizing  function  s0(t,x)  in  (7)  corresponds  to 

the  kernel  of  s0 ,  which  is  given  by 


s0(t,x)  =  Yj^vi(typ,(x)  (31) 

1=1 

Now  note  that  at{t)  =  2yV;(t),  tf>(x)  =  (p(x),  we  see  that  s0(t,  x)  solves  the  optimization  problem 
(1)  since  it  minimizes  the  cost  function  J  ((/>„)  and  at(t) ,  $(x)  satisfy  constraints  (2)  and  (3), 
respectively.  Moreover,  (18)  and  (19)  imply  that  $(x)  is  related  to  at(t)  by 

1  T 

<t>i  (x)  =  —  J  w{t,  x)ai  (i t)dt  (32) 

\  0 

In  the  next  section,  we  show  that  balanced  truncation  is  in  some  sense  similar  to  POD,  in  that,  it 
is  also  optimal  in  the  sense  of  distance  minimization  in  the  Hilbert-Schmidt  norm,  albeit  in 
different  operator  spaces.  The  techniques  developed  for  POD  will  help  us  in  the  context  of 
showing  the  optimality  of  balanced  truncation  as  well. 


ITT.  Optimality  of  Balanced  Truncation 


Balanced  truncation  is  a  simple  and  popular  model  reduction  technique,  which  can  be  described 
as  follows  [21]  [18]:  Suppose  we  have  a  stable  linear  time  invariant  (LTI)  system  described  by 
the  following  ^-dimensional  state  space  equation 

x(t)  =  Ax(t)  +  Bu(t) 

y(t)  =  Cx(t ) 

where  x(t)  is  the  n  x  1  -state  vector  of  the  system,  u(t)  is  an  m  xl -input  vector,  and  y(t)  is  an 
pxl  -output  or  measurement  vector.  A,  B,  and  C  are  constant  matrices  of  appropriate 
dimensions. 


The  underlying  idea  of  balanced  truncation  is  to  take  into  account  both  the  input  and  output 
signals  of  the  system  when  deciding  which  states  to  truncate  with  appropriate  scaling.  The  latter 
is  performed  by  transforming  the  controllability  and  observability  gramians,  denoted  Wc  and 

Wa  respectively,  so  that  they  are  equal  and  diagonal. 


The  controllability  and  observability  gramians  satisfy  the  following  Lypaunov  equations  [18] 

AWC+WCAT  +BBT  =  0  (34) 

AtWo+WoA  +  CtC  =  0  (35) 

The  controllability  and  observability  gramians  can  be  represented  as 

oo  oo 

Wc=  J  eA'BBTeAT,dt,  WQ  =  J  eAT,CTCeA,dt 
0  0 
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Computing  a  state  balancing  transformation  M  is  achieved  by  first  calculating  the  matrix  [18,19], 
Wco  \-WWo ,  and  detennining  its  eigenmodes  Wco  =  M KM  1 . 


z(t)  =  Az(t)  +  Bu(t) 


(36) 


y(t)  =  Cz(t ) 

A:=M~lAM ,  B:=M  lB,  C:=CM  (37) 

The  transformation  M  is  chosen  such  that  the  controllability  and  observability  gramians  for  the 
transfonned  system  satisfy  [18] 

WC=WQ  =  M~xWcM-lT  =  MtWbM  S  (38) 


where  2  is  a  diagonal  matrix  that  satisfies  X2  =  A ,  and  the  diagonal  elements  of  X  ;  cr  ’s,  are 
known  as  the  Hankel  singular  vales  of  the  system,  i.e., 

'L  =  diag{al,a2,  •••,«„}  (39) 


where  cr  are  the  Hankel  singular  values  of  the  system  G  arranged  in  non-increasing  order 

ax  >  a2  >  ■  ■  ■  >  an  >  0  (40) 

In  balanced  truncation  only  states  corresponding  to  large  Hankel  singular  values  are  retained. 
Small  Hankel  singular  values  correspond  to  states  which  are  deemed  weakly  controllable  and 
weakly  observable,  and  therefore  deleted  from  the  state-space  model.  For  instance,  if  the  first  nr 
states  are  retained  then  the  resulting  transformation  is  given  by 

Mr  =  PM  (41) 

where  P  is  the  orthogonal  projection  of  rank  r. 

The  reduced  order  model  is  obtained  by  letting  xr  =  PrMx  as  follows 

xi(t)  =  Aixr(t)  +  Bru(t) 
yr(t)  =  Crxr(t) 

Ar  :=  PrM  ] A MPr;  Br  =PM  ]B,  Cr:=CMPr 
The  error  bound  for  the  output  is  given  by 

\y(t)-yr(t)\2  ^  2Yjai  H0||2 ,  VmgL2  (43) 

nr+ 1 

Balanced  truncation  is  optimal  in  a  precise  sense  when  starting  from  the  so-called  Curtain- 
Glover  balanced  realization.  To  see  this  define  a  causal  bounded  input-output  operator  G  acting 
on  the  standard  space  Z2(- oo,oo)  of  absolutely  square  integrable  functions  defined  on  (-00,00), 
into  Z2  (-oo,  oo)  described  by  the  convolution  [18] 

t 

(< Gu){t )  :=  J  CeA(t  T) Bu{f)dx  (44) 


Now,  define  the  Hankel  operator 
of  G  by 


TG  :  Z2  (-oo,  0]  — »  Z2  [0,  go) 


(45) 
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where  G|z2(_oo0]  denotes  the  restriction  of  G  to  L  (—00  ,0],  and  P+  is  the  orthogonal  projection 

acting  from  L2  (-00,00]  into  Z2[ 0,  00),  i.e.,  P  is  the  truncation  operator 

\f(t)  if  t>{)  , 

o  if  t<0  <46> 

Then,  the  Hankel  operator  TG  can  be  written  as 

0 

r Gu(t)  =  J  CeA(t  T) Bu(j)dz ,  for  t  >  0  (47) 

—00 

The  Hankel  operator  TG  maps  past  inputs  to  future  outputs.  Expression  (45)  shows  that  the 
Hankel  operator  Ify  is  an  integral  operator  mapping  Is  (-00, 0)  into  Is  [0, 00) ,  with  kernel  the 
impulse  response  k(t,r )  defined  by 

k(t,  r)  \=CeA(f~T)B,  r<0,  f>0  (48) 

Balanced  truncation  is  commonly  thought  to  be  a  model  reduction  technique  that  is  not  optimal 
in  any  sense  [22],  We  show  that  this  is  not  the  case,  and  in  fact  balanced  truncation  is  indeed 
optimal  in  the  sense  of  the  Hilbert-Schmidt.  The  techniques  we  use  are  reminiscent  of  the 
previous  section  and  guarantee  for  the  optimum  to  be  a  Hankel  operator.  This  contrasts,  for 
example,  with  the  minimization  in  various  norms  addressed  in  [33]  [34],  To  see  this  note  that  the 
Hankel  operator  Ify  has  finite  rank  k<n  [18],  and  therefore  belongs  to  the  Hilbert-Schmidt 

class  of  operators  acting  from  i}  (—00  ,  0]  into  L  [0, 00)  .  Let  its  spectral  factorization  be  given  by 

Z,sI2(-»,0],  i,€L2l 0,00)  (49) 

1=1 

where  <Ji  are  the  Hankel  singular  values  of  the  system  G  ordered  in  decreasing  order,  i.e., 

al>a2>"->an_1>an  (50) 

and  {x,}]'  and  }[  are  orthononnal  sets  in  f2(-co,0]  and  f2[0,oo),  respectively.  Next, 
consider  the  optimal  distance  minimization 

Pnr  :=  min  1 1  rG  -  r  |fy  (51) 

r  nr<k  r 

where  TG  is  an  operator  acting  from  Z2(-oo,0]  into  f2[0,oo)  of  rank  » .  <  n  .  An  application  of 

identities  (29)  and  (30)  to  the  minimization  (51)  yields  the  unique  optimum  (since  the  distance 
minimization  is  posed  in  a  Hilbert  space) 

r0,=i>,x,®C  (52) 

(=1 

and  the  shortest  distance 

Bn,.  =11  E^-X/  ®  C-  -I>X  ®  c  llasHI  E  °&i  0 C  II HS=  (  E  (53) 

(=1  i= 1  i=nr+ 1  i=nr+ 1 

The  operator  T6  is  not  necessarily  a  Hankel  operator,  however,  we  will  show  that  starting  from 
a  specific  balanced  realization  for  the  original  system,  the  minimizing  operator  can  be  chosen  to 
be  a  Hankel  operator  corresponding  to  the  reduced  order  model.  To  do  so  let  Tc  =  U G{VGT f/1  be 
a  polar  decomposition  of  T0 ,  applying  (16)  and  (17)  to  Tc  the  vectors  x,  and  satisfy 
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(54) 


rGx,  =t/G(nrG)^,=^,.,  i  =  \,-,n 


rGX;  =  UG  (r GrG y  Gi  =  °ili ,  i  =  !,■■■, n  (55) 

That  is,  x,  and  <C  form  a  Schmidt  pair  for  Tc  .  In  terms  of  the  Schmidt  pair  (49)  implies  that  the 
Hankel  operator  Tc  can  be  expressed  as 


(f»(0  = 

i=l 


(56) 


We  propose  the  following  realization  due  to  Curtain  and  Glover  [34]  for  the  impulse  response 
k(t,  t)  given  in  (48),  for  i,  j  =  1, 2,  •  •  • ,  n, 


A  =  (a,):=(^)2f  C;(T)£j(T)dT 

(J  J-00  J 

b  ■-  ( (0),  V^x2  (0),  •  ■  ■  >  x„  (0))r 

C  :=  (^  (0),  (0),  •  •  • ,  (0)) 


(57) 


The  corresponding  semi-group  can  be  computed  as 

a,  -  r0 


eKt=(.-J-Y\  C:(T)Cj(t-T)dr 


(58) 


since 


lim-f°  ^(T)Cj(T-t)dT  =  f  C*(T)C.(T)dT  (59) 

(->0  t  J-00  J  J  — OO  J 

Define  the  controllability  and  observability  operators  denoted  Tc  and  yY 0 ,  respectively  by  [18] 


Note  that  [18] 


x¥c  :Z2(-oo,0]h>M" 

:=  ehT  Bu(t)cIt 
T0  :  M"  i->  Z2[0,oo) 


^oXo’.Ce^x^t^  0 


r  =  t  t 

A  G  A  O  A  C 


(60) 


and  using  the  realization  (57),  we  have 

OF  o'? cu)(t)  =  ZcriCf  £(TMT)dT  =  (TGu)(t) 

J-  oo 

1=1 

and  the  observability  gramian  is  given  by 

*  f 00  *  *  *  ~  7 

=  £  eA  '  C  CeA'<* 

=(J0  (t)Cj(t)dt) 


where  £  is  the  usual  Kronecker  delta. 


/=1  y=l 

=  OA )  =  2  =  diag{a j ,  <r2 ,  •  •  • ,  <7, ) 
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(61) 

(62) 

(63) 

(64) 


Similarly  the  controllability  gramian  TrT*.  =X ,  and  the  realization  (A ,B,C)  is  therefore 
balanced.  By  the  same  token  as  POD  using  a  similar  expression  as  (27),  define  the  Hankel 
operator  corresponding  to  the  nr  -th  order  model  rG  as 

(Tg  u){t)  =  r)dr 

(65) 

=  f  C  PreM'~r) Pr  B  dr 

J  —00 

(66) 

=  f°  C  Pr  {PreMt~z)Pr  )PrBdr 

J— 00 

(67) 

The  last  equality  follows  by  (57),  (58)  and  the  fact  that  P2  =  Pr  (since  Pr  is  a  projection).  Putting 


Cr  :=CPr,  Br:=PrB  correspond  to  truncating  C  and  B,  respectively,  and  (59)  implies  that 
PeM'-:] Pr  =  ,  and  Ar  :=  Pr  A  P  correspond  to  truncating  the  state  space  model 

(A ,B,C)  to  nr  states,  and  the  Hankel  operator  has  rank  nr .  Moreover, 

ftir=l|r0-r(,J|„=(  £  f  (68) 

i—nr  +1 

By  uniqueness  of  the  minimizer  in  (51),  expressions  (65)  and  (68)  imply  that  we  must  have 

r  =  r 

L  r  -  1  G„r  • 

In  terms  of  kernel  approximation,  balanced  truncation  is  a  particular  case  of  POD  in  the  sense 
that  the  kernel  we  want  to  approximate  is  the  impulse  response  of  the  system  k(t,  r)  defined  in 

(51).  The  optimization  index  /un  can  then  be  written  as  in  POD 

Pn,  =min{£  f °Jk(t,T)-YJfi(t)gi(T)\2  drdP.f,  eT2[0,co);g.  eT2(-oo,0]}  (69) 

i=l 

=  r  [°  I  k(t,  r)  -  C  eAr(‘~T>  Br  I2  dzdt  (70) 

JO  J— oo 

Expressions  (54)  and  (70)  show  that  balanced  truncation  is  optimal  in  the  sense  of  optimal 
approximation  in  the  Hilbert-Schmidt  norm  of  the  Hankel  operator  TG ,  and  optimal  in  the  sense 

of  the  IHIj-norm  of  kernels  corresponding  to  impulse  responses  of  linear  time-invariant  systems 
defined  over  [0,co)x(-cx>,0] .  The  linear  time-invariant  system  framework  allows  the  exact 
computations  of  the  optimal  lower  order  model  approximation.  This  contrasts  with  POD  which 
uses  simulation  data  and  particular  open-loop  inputs  to  generate  snapshots. 

IV,  Prototype  Problem:  Nonlinear  Convection 

The  specific  problem  geometry  considered  is  shown  in  Figure  1.  The  idea  and  methods  presented 
could  be  modified  to  apply  to  a  different  geometry  or  obstacle  shape.  A  realistic  example  of  this 
geometry  in  an  aerodynamic  application  is  a  payload  hatch  open  during  flight  with  actuator 
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control  only  on  the  boundary.  Let  Q^be  the  region  defined  by  [ai,  a2]  x  [bi,  b2].  Let  ClfllU  be  the 
region  defined  by  (a0,  aend)  x  (bo,  bend).  Then  the  problem  domain  is  given  by  Q  =  Clfull  /  Q  .  In 
this  problem  setup,  Qgapis  an  obstacle.  The  system  dynamics  that  act  within  the  problem  domain 
are  described  by  the  two  dimensional  (2D)  Burger’s  equation 


d  1 

—  w(t,  x,  y )  +  V«F(w)  =  - 
dt  r 


^2  ^2 

w(t,x,y)  +  —rw(t,x,y) 


dx 


where  the  form  of  F(w )  is 


F(w)  = 


w\t,x,y) 


dy 2 


w2(t,x,y) 


(71) 


(72) 


The  value  r  is  similar  to  the  Reynolds  number  used  in  the  Navier-Stokes  Equation.  This 
parameter  controls  how  much  nonlinearity  is  present  in  the  problem.  The  2D  Burgers’  equation 
is  chosen  as  a  surrogate  for  the  Navier-Stokes  equations,  since  it  has  a  similar  convective 
nonlinearity  but  can  be  coded  on  a  PC. 


Figure  1.  Problem  Geometry 


The  boundary  conditions  on  the  top  and  bottom  are  described  by  the  following  equations 

W(/,r  bottom)  =  UbottorSf)^  bottom  (*)  (73) 

^t,T,op)  =  utopmtop(x)  (74) 

Here  w,0/J(t)and  uhottom(t)  are  the  control  inputs  on  the  top  and  bottom  boundaries  respectively, 
the  spatial  functions  4f  (x)  and  vP6o(tom(x)  describe  the  spatial  effect  that  the  controls  have  on 
the  top  and  bottom  boundaries.  The  boundary  condition  on  the  airflow  intake  side  is 
w(t,  r  )  =  f(y)  and  is  parabolic  in  nature.  The  airflow  outtake  side  has  a  Neumann  boundary 

d 

condition  that  has  the  form  — w(t,T ' out)  =  0  .  On  all  of  the  remaining  boundaries  of  Cl ,  w(t,x,y ) 

dx 

is  set  equal  to  0  for  all  values  of  t.  Finally,  the  initial  conditions  for  the  interior  are  given  by 
w(0,x,y)  =  w0(x,y )  e  L2(Cl). 
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A  numerical  solution  was  found  by  simulation  using  a  uniformly  spaced  grid.  The  resulting 
system  model  contains  a  little  more  than  2000  states.  After  this  a  method  of  optimal  (with  respect 
to  average  kinetic  energy  [27])  subspace  detection  was  performed.  POD  is  a  method  for  model 
reduction  that  works  on  the  idea  that  a  set  of  basis  vectors  in  an  infinite  dimensional  space  can  be 
created  that  gives  the  best  representation  of  typical  system  behavior  for  a  fixed  basis  dimension 
size  [27].  The  method  used  for  determining  how  many  modes  need  to  be  retained  is  a  design 
choice.  A  reasonable  choice  and  the  one  used  in  this  POD  model  construction  were  to  look  at  the 
total  energy  captured.  When  a  vector  having  significant  system  “energy”  was  found,  it  was 
forced  to  be  perpendicular  to  all  previously  found  <pk' s .  A  condition  that  99.9%  of  system 

energy  must  be  captured  was  used  for  detennining  how  many  system  modes  were  needed.  This 
condition  was  met  by  a  40  POD  basis.  Although  this  is  a  major  reduction  from  the  numeric 
solution,  it  will  be  shown  that  important  system  dynamics  can  be  retained  with  even  lower  state 
number  system  models. 

Other  names  for  this  decomposition  are  the  Karhunen-Loeve  decomposition  [36][37],  and 
principal  component  analysis  [37].  Some  applications  have  been  in  fluid  mechanics,  random 
variables,  image  processing,  and  data  compression  [10].  This  decomposition  method  has  been 
used  extensively  ([1]  [4][8][7][11][12][13][14]  [15][38][39][40][41])  to  analyze  experimental  or 
simulation  data  and  to  extract  the  most  dominant  trends.  The  general  approach  of  this  method  is 
to  construct  a  series  of  solution  “snapshots.”  These  snapshots  are  solved  by  numeric  simulations 
of  the  governing  system  equation(s)  with  a  variety  of  input  equations.  The  method  of  finding 
these  solution  snapshots  is  described  in  detail  in  [42]  [43].  The  greatest  strength  of  this  method  is 
that  it  provides  a  relevant  and  optimal  set  of  basis  functions  that  allows  a  low-dimensional 
subspace  to  be  identified.  This  subspace  allows  us  to  derive  a  model  when  the  governing 
equations  are  projected  onto  the  subspace  [27].  The  snapshots  are  needed  to  generate  the 
correlation  matrix  that  can  be  used  to  find  the  POD  basis.  It  is  important  to  choose  relevant  input 
signals  for  the  numerically  simulated  system.  Further,  these  inputs  should  be  similar  to  the 
expected  inputs  of  the  real  system.  In  this  way,  the  POD  model  should  accurately  represent  the 
system  in  the  nonnal  operating  range.  In  this  paper  the  inputs  used  for  system  identification  are 
of  the  form  [12] 


ubot,om(  0  =  ^sin(0.25t2)  ulop(t)  =  0,  (75) 

ubot,om(  0  =  0  *V(0  =  /?sin(0.25t2),  (76) 

ubottom(t)  =  /?sin(0.25r)  utop(t)  =  /?sin(0.25 12)  ,  (77) 

where  the  values  for  J3  are  -3,  -2,  and  -1  and  the  range  for  t  is  0  to  10  seconds  with  a  sample 
every  50  milliseconds.  The  squelch  signal  for  all  three  values  of  f5  is  shown  in  Figure  2. 

The  numerical  simulation  was  performed  to  create  the  ensemble  of  solution  snapshots 
[Sk(x,y)}Mk=i  [12].  The  snapshots  refer  to  a  collection  of  samples  at  specific  points  of  individual 
solutions  to  the  governing  equation  [27]. 

Each  snapshot  captures  samples  at  specific  points  of  the  problem  geometry  while  the  control 
inputs  are  varied.  The  value  for  M  (the  number  of  snapshots)  must  be  greater  than  the  number  of 
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modes  that  one  will  choose  for  the  approximated  system  model.  For  a  good  representation  M 
should  be  much  larger  than  the  desired  size  for  the  POD  basis  [12]. 


Inputs  for  System  Excitation 


Figure  2.  Test  Inputs  Used  to  Generate  for  Snapshots 


The  solution  to  the  PDE  is  assumed  to  be  finite  energy,  i.e.,  belongs  to  the  Hilbert  space 
L  ([(),  7’)xQ)  .  The  solution  is  approximated  as 


w(t,x,y)*s'£ak(t)0k(x,y)  where  akeL\[0,T)),  </>k(x,y)eL2(Cl)  (78) 

k=\ 

where ak  (t)  eL2{[0,T))  are  time  varying  coefficients  that  multiply  the  spatial  functions 

</>k(x,y)eL2( Q),  where  T2([0,  T))  and  L2(fl)  are  the  standard  Hilbert  spaces  of  absolutely 
square  integrable  functions  defined,  respectively,  on  the  time  interval  [0,T)  and  spatial  domain 
Q  .  The  approximation  (78)  can  be  as  accurate  as  desired  since  the  tensor  space 

L2([0,T))  ®L2 (Q)  :=  | Yjxk(t)(/)k{x,y),  ak(t)  eL2([0,T)),  <f>k(x,y)  gT2(Q),  V»  integerj  (79) 


is  dense  in 


L  ([0,7>O)  .  Any  basis  for  L2(fl) 


can  be  used  to  construct  the  approximation  of 


the  solution  w(t,x, y) .  Here  we  use  the  POD  basis  {<f>k}  since  it  was  shown  to  be  optimal  in  the 
following  sense 


Mn  :=  inf 

ak^L\[0J)\ 

4(x,y)ei2(U) 


w{t,  x ,  y )  -  Yj  ak  (0  A  (x,  y) 


k= 1 


and 


A, 


->0 


(80) 


where  ||  •  ||7  is  the  norm  of  L  ([0, 7)xfij .  In  section  II  it  was  shown  that  (80)  corresponds  to 


optimal  approximation  of  an  integral  operator  with  kernel  w(t,x,y)  in  the  HS  nonn  by  finite 
rank  HS  operators.  A  numerical  solution  was  found  by  simulation  using  a  uniformly  spaced 
grid.  The  resulting  system  model  contains  a  little  more  than  2000  states.  The  POD  model 
construction  is  based  on  the  total  energy  captured.  A  condition  that  99.9%  of  system  energy 
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must  be  captured  was  used  for  determining  how  many  system  modes  were  retained.  Although 
this  is  a  good  reduction  from  the  numerical  solution,  it  is  shown  that  important  system  dynamics 
can  be  retained  with  even  lower  order  models.  The  general  approach  of  POD  is  to  construct  a 
series  of  solution  “snapshots.”  These  snapshots  were  generated  by  numerical  simulations  of  the 
2D  Burgers’  equation  with  a  variety  of  chirp  input  signals.  Construction  of  a  POD  basis  capable 
of  spanning  the  baseline  solution,  as  well  as  dynamics  introduced  by  boundary  actuation  has 
been  addressed  using  split  POD.  A  Galerkin  type  projection  using  Divergence/Green's  theorem 
results  in  a  nonlinear  Galerkin  model.  The  Galerkin  projection  provides  only  a  weak  solution. 
However,  this  weak  solution  with  finite  difference  approximations  of  the  boundary  conditions 
eventually  leads  to  a  nonlinear  temporal  model  for  the  temporal  or  POD  coefficients.  This 
solution  does  not  explicitly  include  the  control  inputs  or  boundary  condition  information  into  the 
governing  equation.  In  order  to  do  so  an  approximation  of  the  partial  derivatives  is  carried  out 
including  the  control  inputs  and  the  boundary  data.  After  substitutions  w(t,  x,  y )  can  be 

approximated  as  a  linear  combination  of  POD  modes  when  the  ak(t)'s  are  solved  in  the 
nonlinear  Galerkin  model. 

The  ultimate  goal  for  finding  a  solution  to  the  governing  equations  is  then  to  find  the 
aks  and  <pk's .  The  correlation  matrix  L,  of  size  M  x  M,  is  composed  of  the  inner  products  of  the 
snapshots.  Then  in  L 2  (Q) 

LtJ=(S„S,)  (81) 

(S.^^jj^dxdy  (82) 


where  *  denotes  complex  conjugate  transpose. 

The  term  ‘ensemble  member’  is  used  to  describe  one  element  of  the  collection  of  solution 
snapshots  [12]. 


The  function  (pt  is  chosen  to  maximize  the  average  projection  of  the  member  ut  onto  q>r  This 
idea  is  expressed  in  the  following  manner  [27]  [28] 

Wi)  r 


max 

q?<=L2  (Q) 


(83) 


A  choice  must  be  made  on  how  many  modes  will  be  kept  in  the  POD  constructed  subspace.  This 
number  can  be  chosen  to  be  sufficiently  large  such  that  the  majority  of  the  system’s  total  energy 
is  captured.  Additionally  this  number  of  modes  to  keep  could  be  based  upon  the  hardware 
abilities  of  the  device  doing  the  computations.  The  chosen  number  of  modes  to  keep  is  denoted 
as  n.  A  singular  value  decomposition  of  the  matrix  L  is  performed.  The  n  largest  eigenvalues 
[\,A2,.. ,,An)  of  the  matrix  L  are  found  and  placed  in  descending  order.  Then  the  set  of 

eigenvectors  are  identified  to  be 


16 


The  resulting  orthonormal  POD  basis  of  dimension  n  can  be  constructed  using  the  information 
found  from  the  correlation  matrix  L.  First,  the  eigenvectors  of  L  are  weighted  by  their 
corresponding  eigenvalues  and  normalized  according  to  [12] 

K  ||VA- 1|2  =  1 ,  for  k  =  {l,2,...,«}  .  (84) 

Then,  the  POD  basis  set  is  formed  according  to 


M 

<pk(x,y)  =  £vkJsj(x,y X  (85) 

j= i 

with  vk  .  being  the  jth  component  of  the  eigenvector  vk  .  Solving  equation  (85)  gives  the  n  <pk's 
which  constitute  the  POD  basis  of  dimension  n. 


The  governing  equation  is  projected  onto  the  POD  basis.  The  projection  is  accomplished  via  a 
Galerkin  type  projection.  The  Galerkin  projection  results  in  only  a  weak  solution  to  the  PDE. 
However,  this  weak  solution  with  finite  difference  approximations  of  the  boundary  conditions 
eventually  leads  to  a  nonlinear  temporal  model  for  the  temporal  or  POD  coefficients  {ak}  [4], 
Projecting  (1)  onto  the  POD  basis  yields  [4] 


0  1 

J  —w{t,x,y)(pk{x,y)dx  =  —  ( J  (Vw(t,x,y)*n)^(x,y)^(x)-| Ww{t,x,y)>W  (pk{x,y)dx)- 

Q  ^  dCl  Cl 

( |  (F(w)‘n)<pk  (x,  y)dA(x)  -  F-V  (pk  (x,  y)dx) 


dCl 


where  the  first  term  on  the  right  hand  side  is 

d 


f  (Vw(f  x,  y)»n)(pk  (x,  y)dA{x)  =  f  —  w(t,  x,  bx  )<pk  (x,  b]  )dx  -  f  ‘  —  w(t,  x,  b2  )<pk  (x,  h,  ))dx 

J  J  Pi?  '  Ja , 


dCl 


8y 


r 4~w(t,  a0 ,  y)<pk  (a0 ,  y)dy 
Jh°  ox 

The  Neumann  boundary  condition  forces  the  portion  of  the  boundary  integral  over  bo  to  bend 
along  aend  to  be  0,  i.e., 

f""  fMt,  “end  >  y)(Pk  (ae„d  »  =  0  . 

Jh0  OX 

The  second  boundary  integral  is  decomposed  as  follows 

j*  1  j*  ^ snd  2  2 

I  , (F(w)»n)<pk (x, y)dA(x)  =  -  (w(t,aend, y)  <Pk(aend, y)-f(y)  <pk(a0,y))dy 

This  solution  does  not  explicitly  include  the  control  inputs  or  boundary  condition  information 
into  the  governing  equation.  In  order  to  do  so  an  approximation  of  the  partial  derivatives  is 
carried  out  including  the  control  inputs  and  the  boundary  data.  If  h  denotes  the  step  size  between 
the  points  on  the  unifonn  Cartesian  grid  used  for  the  finite-difference  solution,  then  we  have  [12] 
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x  b ) ~  -/?) 

dy  ’  ’  1  h 


d  ,  ,  ,  w(t ,x,b2+h)~  u  (t )VF  (x) 

—  w(f  x,b2)~ - - - 

dy  h 


d 

—  w(t,a0,y ) 
ox 


w(t,a0+h,y)-f(h ) 
h 


After  substitutions  w(t,x,y)  can  be  approximated  as  a  linear  combination  of  POD  modes  when 
the  ak '  s  are  solved  in  the  following  system  model.  Then,  the  temporal  model  for  the  system  is 
given  by  [12] 

a  =  Aa  +  Bu  +  N(a)  +  F,  a(0)  =  ao  (86) 

where  a  e  M"  and  the  matrices  A  is  n  x  n,  B  is  n  x  2,  N  and  F  are  both  vectors  n  x  1 . 

The  output  equation  will  be  simply  chosen  to  be 

y(t)=cc(t) 

In  this  model  the  dimension  of  the  state  vector  a  is  40  which  correspond  to  40  POD  modes.  The 
first  8  POD  modes  corresponding  to  the  first  8  temporal  coefficients  are  shown  in  Figure  3.  The 
first  model  corresponds  to  the  baseline  mode,  and  the  remaining  modes  to  actuated  modes. 


In  Figure  4,  dashed  lines  denote  the  linear  combination  of  POD  modes  restricted  to  the  boundary. 
Solid  lines  denote  the  boundary  test  inputs  defined  by 


ux(t)= sin 


2nt 


w2(t)  =  sin 


2nt 


As  can  be  seen  in  Figure  4,  there  is  very  good  agreement  between  the  boundary  conditions 
specified  for  the  full  order  system  and  the  linear  combination  of  POD  modes  restricted  to  the 
boundary. 

The  goal  of  model  reduction  is  to  construct  another  nonlinear  system  [1 1  ]  [5] 


ar  =  Aar  +  Bru  +  Nr  (ar)  +  Fr 


(87) 


where  ar  e  Mr  and  r  <  n  ,  such  that  the  behavior  of  the  two  systems  is  similar  for  states  in 

some  region  of  the  state  space.  The  reduced  model  is  derived  via  the  construction  of  an 
immersion/projection  pair 


a  =  far  ar=Ta  TT  =  Ir 


(88) 


where  I r  is  the  rxr  identity  matrix,  resulting  in  the  following  reduced  model 
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Figure  3:  First  Eight  POD  modes 


The  top  control 


The  bottom  control 


Figure  4:  Boundary  Control  Accuracy 


19 


ar  =  TATar  +  TBu+  TN(Tar)+  TF 

y(t)  =  CTar(t ) 


(89) 


This  is  carried  out  by  developing  an  empirical  balanced  truncation  algorithm  which  is  based  on 
experimental/simulation  input-output  measurements  of  the  nonlinear  Galerkin  model.  This  is 
introduced  in  the  next  section. 


V.  A  New  Empirical  Balanced  Truncation 

Balanced  model  reduction  requires  the  knowledge  of  the  controllability  and  observability 
gramians.  The  latter  are  obtained  by  solving  Lyapunov  equations,  which  is  prohibitive  for  large 
scale  systems.  For  a  system  with  n  states,  the  controllability  and  observability  matrices  are 
symmetric  matrices  and  therefore  solving  each  one  of  them  involve  finding  n(n+l)/2  unknowns. 
An  alternative  is  to  develop  a  balanced  truncation  algorithm  based  on  empirical  gramians,  which 
are  constructed  solely  from  a  single  simulation  using  a  sufficiently  rich  input.  The  proposed 
empirical  balanced  truncation  uses  the  Galerkin  model  with  chirp  signals  as  inputs  to  produce  the 
outputs  in  the  so-called  Eigensystem  Realization  Algorithm  (ERA).  The  latter  estimates  the 
system’s  Markov  parameters  that  accurately  reproduce  the  outputs.  Balanced  truncation  is  used 
to  show  that  model  reduction  is  still  effective  on  ERA  produced  approximate  systems.  A  new 
method  of  finding  empirical  controllability  and  observability  gramians  for  the  approximated 
system  is  introduced.  After  the  empirical  gramians  are  approximately  balanced,  the  necessary 
transformation  maps  are  applied  back  to  the  original  system.  The  empirically  balanced 
realization  is  then  truncated  to  further  reduce  the  system  model  size  while  retaining  the  most 
important  system  dynamics.  The  linear  empirical  balanced  truncation  algorithm  is  applied  to  the 
Galerkin  model  which  is  nonlinear.  The  rationale  for  doing  so  is  that  linear  subspace 
approximations  to  exact  submanifolds  associated  with  nonlinear  controllability  and  observability 
require  only  standard  matrix  manipulations  utilizing  simulation/experimental  data. 


The  computation  cost  to  solve  large  Lyapunov  equations  for  the  controllability  and  obervability 
gramians  prompt  us  to  propose  a  balanced  truncation  algorithm  based  on  empirical  gramians 
constructed  from  input-output  data  measurements.  To  this  end,  let  us  first  introduce  the  /-step 
observability  and  g-step  controllability  matrices  [47] 


o,-= 


c 

CA 

CA2 


,  Rq  :=[B  AB  A2B  •••  AqAB~\ 


(90) 


CA'-1 


which  give  rise  to  the  /-step  observability  and  g-step  controllability  gramians 


w  =0*0 

''of 


wcq-=RqRl 


(91) 


As  the  numbers  q  and  /  approach  infinity,  these  empirical  gramians  approach  the  true  gramians 
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lim  Wol=Wa,  lim  Wcq=We  (92) 

The  goal  is  to  find  a  balancing  transformation  matrix  M  that  will  approximately  balance  the 
empirical  gramians,  i.e., 


Wcq  :=  MWcqM*  =  (M*yxWolM~l  =:  Wol  =  X  (93) 

The  matrix  M  can  then  be  applied  back  to  the  original  system  model  to  produce  an  approximately 
balanced  realization. 

The  product  of  the  /-step  controllability  and  the  c/-step  observability  matrices  gives  a  Hankel1 
matrix,  denoted  Hlq ,  containing  the  Markov  parameters  C/1'  B,  k  =  0,1,  •••  ,  of  the  system  in  the 
following  way 


H,q  :=  0,R  = 


CB 

CAB  • 

■■  CAq  B 

CAB 

CA2B  ■  ■ 

■■  CAqB 

CAllB 

CA'B  ■  ■ 

■■  CA1+qlB 

(94) 


Y/>  1 


1 

M 

o 

1 

* 

1 _ 

1 

O 

- 1 

o 

_ 1 

for  integers  /  and  q  chosen  such  that  [47] 

rank(Hlq )  =  rank(H (l+m+j)  )=n. 

In  terms  of  the  SVD  decomposition  of  HIq 

H,q  =cxr  =[C,  U2] 

The  balancing  transformation  M  is  constructed  as 

m=r,v,  r/> 

A  straightforward  computation  shows 

Wcq  :=  M-lWcqM*~'  =  M*WolM  =:  Wol  =  X, 

Balanced  truncation  can  be  realized  the  usual  way,  if  cr  »  <Jr+x  for  some  r  then  we  can  partition 
Xj  as 

X  o 

o  X,. 


(95) 

(96) 

(97) 

(98) 


2:= 


r+\ . 


(99) 


where 


1  A  Hankel  matrix  is  simply  a  matrix  that  has  the  ih  column  identical  to  the  ith  row. 
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(100) 


Z,  =  diag(c Tva2,---,ar),  Er+1  =  diag(ar+l,ar+2,---,an) 

A  columwise  confonnal  partition  of  Ux  and  Vx 

Uy=\Ur  Un_r],  Vx=[Vr  Vn_r ]  (101) 

yields  the  immersion/projection  pair 

f=RV'L~^2,  T  =  Y  ^  U*  O,,  TT=I  (102) 

and  from  which  a  reduced  order  r-dimensional  model  with  state  matrices  is  deduced 

Ar  =TAT  Br  =  TrB  Cr=CTr  (103) 

The  above  construction  only  requires  estimates  of  the  Markov  parameters 

CAkB,  k  =  0,\,---,l  +  q-l  (104) 

A  basic  relationship  between  the  Markov  parameters  and  the  input  and  output  relationship  in 
discrete-time  is 

oo 

y(k)  =  XT(f>w(k-f>  (105) 

(=0 

7(0)  =  D,  7(1)  =  CB,  7(2)  =  CAB,  •  •  • ,  Y(k)  =  CAk~lB  (106) 

The  Markov  parameters  can  be  computed  from  a  single  simulation/experiment  in  which  a 
sufficiently  rich  input  signal  is  applied  and  the  output  responses  are  collected.  In  the  next  section, 
the  Discrete  Fourier  Transform  (DFT)  is  used  to  map  time  domain  data  into  spectral  densities 
from  which  frequency  response  estimates  are  calculated  using  the  Eigensystem  Realization 
Algorithm  (ERA)  [44], 

VI,  Eigensystem  Realization  Algorithm  (ERA) 

Several  frequency  domain  identification  techniques  are  used  in  practice  to  identify  the  model 
parameters.  One  such  method  that  is  a  well-defined  frequency  domain  system  identification  is 
the  Eigensystem  Realization  Algorithm  (ERA)  technique  [44].  The  ERA  based  system 
realization  model  is  created  directly  from  empirical  data  and  frequency  domain  characteristics  of 
transfer  functions.  This  method  is  applied  to  discrete  time  versions  of  system  models. 

The  ERA  and  the  empirical  balancing  method  have  been  in  many  cases  [48][49][50][5]  to 
identify  the  system  based  upon  experimental  or  simulated  results.  The  workhorse  of  the  method 
involves  finding  the  system  Markov  parameters  to  populate  a  matrix  that  gives  the  impulse 
response  from  each  input  to  each  output  at  various  time  steps. 

The  ERA  process  starts  by  forming  the  pulse  response  matrices  Y(k).  Theoretically,  pulse  inputs 
are  sufficient  to  excite  the  system  to  find  the  approximate  frequency  domain  transfer  function. 
The  impulses  must  be  timed  so  that  they  do  not  all  occur  at  the  same  time  so  the  individual 
impulses  can  overlap  and  are  distinct  for  each  input  channel.  However,  practically  a  single 
impulse  does  not  have  enough  energy  to  fully  excite  every  state  [50].  Each  additional  pulse 
response  matrix  is  calculated  such  that  Yt  j{k)  is  the  ith  output  at  sample  k  resulting  from  a  pulse 

input  on  the/7'  input  channel.  The  dimension  of  each  of  the  Markov  parameter  matrices  is  p  x  m 
where  p  is  the  number  of  system  inputs  and  m  is  the  number  of  system  outputs.  As  k  increases, 
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the  size  of  the  Hankel  matrix  composed  of  the  Markov  parameters  increases.  This  matrix  does 
not  have  to  be  square.  In  practice,  a  broadband  random  input  or  a  chirp  is  applied  and  estimates 
of  the  input  autocorrelation  and  input/output  cross-correlation  are  used  to  identify  the  system 
Markov  parameters  [44], 

An  alternative  form  to  (105)  can  be  created  not  using  the  actual  outputs  and  inputs  but  replacing 
the  output  term  by  the  cross-correlation  between  the  inputs  and  the  corresponding  outputs  [44] 

oo 

=  (107) 

8=o 

where  the  length  of  the  data  sequence  is  m, 

1  m- 1  1  m- 1 

Ru,xk)=—Tju^uT^k~r^  and  Rvu(k)=—YJy(^uT(k~T^  (i°8) 

m  z--o  '  m  r~Q 

The  basic  process  for  finding  the  Markov  parameters  starts  using  the  ratio  of  the  power  spectral 
density  of  the  cross-correlation  between  the  inputs  and  outputs  and  the  power  spectral  density  of 
the  autocorrelation  between  the  input  signals. 

These  power  spectral  densities  are  given  by  the  following 


i  m- 1  _  ;2?rfc 

pvu(k)=— IX»  m  and 

m  r=0 


1 

W  =  -^Rm(r)e 

m  r=0 


.Ink 


(109) 


The  ratio  of  the  two  power  spectral  densities  is  the  frequency  response  function  and  is  given  by 
[44] 

P  (AT 

^tiry  <110) 

Then,  the  final  step  is  to  take  the  inverse  Fourier  transform  to  find  the  pulse  response  (Markov 
parameter)  matrices  [44] 


°°  ] 

Yk:=  Y(k)  =  ^G(zk)e 

£=0 


(Hi) 


This  system  identification  method  can  give  good  results  even  in  the  case  where  very  noisy  data  is 
acquired.  This  ability  lies  in  taking  many  more  Markov  parameters  than  the  true  system  order. 
The  computed  Markov  parameters  that  do  not  agree  with  the  other  data  can  effectively  be  filtered 
out  of  the  frequency  response  function.  The  Hankel  matrix  containing  the  Markov  parameters  is 
of  the  following  form  [67] 


H, 


Y  Y 
1\  1 2 

Y  Y 
12  A 


F 


F 


<7+ 1 


(112) 


F  F 


/+1 


V 


q+l-l 


The  individual  Ft’s  correspond  to  the  following  sequence: 
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Y0  =D,  Yx=  CB,  Y2  =  CAB,  Yk  =  CAk  lB 


(113) 


In  some  cases  the  input  data  for  the  ERA  method  might  be  provided  by  an  experiment  on  a  real 
system.  However,  in  this  project  a  unique  approach  of  using  the  Galerkin  model  in  the  place  of 
the  real  system  was  used  to  generate  the  empirical  data.  The  full  order  system  model  was  created 
using  finite-difference  methods.  Recall  that  the  control  inputs  were  explicitly  placed  in  the 
boundary  conditions  because  the  control  inputs  do  not  show  up  explicitly  in  the  two-dimensional 
burger’s  equation.  However,  the  weak  Galerkin  model  results  in  a  nonlinear  state  space  model 
that  simplifies  the  relationship  between  the  input  and  outputs. 

The  chirp  signals  used  for  the  excitation  of  the  Galerkin  model  are  of  the  following  form  and  are 
shown  in  Figure  5. 

w,(t)= -sin(0.55t2)  (114) 

u2(t)  =  -sin(0.60t2)  (115) 


ERA  Method  Excitation  Control  Inputs 


Figure  5.  Excitation  Inputs  for  ERA  Method 


VII,  Application  to  the  Galerkin  Model 

The  empirical  balanced  truncation  based  on  linear  systems  is  applied  to  the  Galerkin  model 

a  =  Aa  +  Bu  +  N(a)  +  F,  a(0)  =  ao  (116) 

which  has  an  equilibrium  in  steady  state,  denoted  by  ass .  The  rationale  for  doing  so  is  that  linear 

subspace  approximations  to  exact  submanifolds  associated  with  nonlinear  controllability  and 
observability  require  only  standard  matrix  manipulations  utilizing  simulation/experimental  data 
as  explained  in  [49].  The  computational  advantages  of  the  scheme  presented  here  carry  over 
directly  to  the  nonlinear  setting. 
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The  reduced  order  model  is  derived  as  discussed  through  the  construction  of  the 
immersion/projection  nonlinear  system  pair 


a=Tar,  ar  =  Ta 


(117) 


This  results  in  the  following  reduced-order  model 


a,  =  Arar  +  Bru  +  Nr  (T  ar  )  +  Fr,  ar  (0)  =  Ta0 


(118) 


where 


4  :=  TAT, 


Br  :=TB,  Nr  :=TN,  Fr  :=TF 


If  (1 16)  has  a  linearization  around  the  steady  state  equilibrium  ass 

a  =  A(a  +  Beu,  a(0)  =  ao 


where 


A.  = 


4  = 


d(Aa  +  Bu(t)  +  N(a)  +  F) 


da 


d(Aa  +  Bu(i)  +  N(a)  +  F) 


du 


and  the  reduced  system  linearization  around  the  steady  state  equilibrium  as! 

ccr  =  Atrar  +  B(ru,  ar  (0)  =  T  ao 


where 


4,= 


d(Arar  +  Bru(t )  +  Nr  (f  ar )  +  Fr ) 


4,  = ' 


da 


d(Arar  +  Bru{t)  +  Nr{Tar)  +  Fr) 


du 


(119) 

(120) 

(121) 

(122) 

(123) 

(124) 


The  linearization  of  both  models  about  equilibrium  ass  are  related  by 

Ar  =  TAtf,  B(r  =TB( 

Empirical  balanced  truncation  applied  to  the  40  order  Galerkin  model  resulted  in  a  14th  order 
reduced  model.  The  first  8  temporal  coefficients  of  the  14  order  reduced  model  and  2000  full 
order  model  are  plotted  in  Figure  6.  Figure  6  shows  good  agreement  between  the  temporal 
coefficients. 


In  Figure  7,  we  compare  the  Hankel  singular  values  of  the  2000  full  order  linearized  and  reduced 
14th  order  empirical  model.  As  expected  the  Hankel  singular  values  corresponding  to  the 
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reduced  order  model  are  smaller  than  the  full  order  model,  nevertheless  the  Figure  shows  that 
they  are  close.  In  Figure  8  we  compare  the  full  order  solution  w(t,x,y)  of  the  Burgers’  equation 

with  the  solution  based  on  the  14th  order  ERA  model  vr.(i,x,_y)  The  Figure  shows  that  they 
behave  similarly  especially  at  the  boundary  where  control  is  applied. 


Non  Linear  Galerkin  Model  and  ERA  Model  Truncated  to  14  States 
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Figure  6.  Projected  and  POD  Model  Coefficients. 
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Figure  7.  Comparison  of  Flankel  Singular  Values. 
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Full  Solution 
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Reduced  Solution 


Figure  8.  Full  and  Reduced  Order  Models’  Responses. 

VIII.  A  New  Empirical  Hankel  Norm  Model  Reduction 

In  this  section,  we  develop  and  apply  an  empirical  Hankel  norm  model  reduction  technique  to  a 
Galerkin  model.  The  latter  is  obtained  by  projecting  the  governing  equations  of  a  nonlinear 
convective  flow  on  its  proper  orthogonal  decomposition  basis.  The  nonlinear  convective  flow  is 
chosen  as  a  prototype  problem. 

The  proposed  empirical  Hankel  norm  model  reduction  uses  the  Galerkin  model  with  a  chirp 
signal  as  input  to  produce  an  approximate  linear  model  and  estimates  its  Markov  parameters  that 
accurately  reproduce  the  output.  A  method  of  finding  empirical  controllability  and  observability 
gramians  for  the  approximated  system  is  discussed.  After  the  empirical  gramians  are 
approximately  balanced,  the  necessary  balancing  transformation  matrix  can  be  applied  back  to 
the  original  system  to  get  an  approximate  balanced  realization.  Further,  Hankel  optimal  norm 
reduction  is  applied  to  the  latter  model.  Hankel  model  reduction  uses  the  system  operator 
induced  norm  instead  of  the  average  system  energies  used  in  POD  and  balanced  model  reduction 
techniques.  The  operator  norm  is  a  more  precise  and  appropriate  measure  of  system  modeling 
errors.  Application  to  the  nonlinear  convective  flow  as  a  prototype  problem  shows  the 
effectiveness  of  the  proposed  empirical  Hankel  nonn  model  reduction,  and  promise  for 
application  to  the  Navier-Stokes  equations  since  the  fonner  has  a  similar  nonlinearity. 
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Hankel  norm  approximation  is  an  optimal  model  reduction  in  the  sense  of  the  operator  norm  of 
the  Hankel  operator  associated  to  the  system  [51][18][19][52][53].  This  makes  it  particularly 
useful  for  model  reduction  since  it  is  desirable  for  the  error  induced  noun  between  the  original 
and  reduced  model  be  small.  The  idea  is  to  approximate  the  Hankel  operator  associated  to  the 
system  by  another  Hankel  operator  with  lower  rank  in  an  optimal  way.  The  error  induced  norm  is 
the  worst  case  error  due  to  the  input  and  making  it  small  guarantees  that  the  input-output 
behaviors  of  both  the  original  and  reduced  models  are  close  for  any  input.  In  practice  the 
algorithm  uses  a  balanced  state  space  representation  [18]. 

Hankel  norm  approximation  provides  a  better  error  bound  than  balanced  truncation  and  works 
well  even  when  there  is  no  Hankel  singular  value.  However,  there  is  no  viable  Hankel  norm 
model  reduction  for  nonlinear  systems  and  we  resort  to  constructing  an  approximate  linear  model 
from  empirical  data  using  the  nonlinear  Galerkin  model. 


Recall 


(Gu)(t)  *=  f  CeA(t~T) Bu(t)<It  =  y{t) 

J  —oo 


(125) 


The  Hankel  operator  rG  can  be  written  as 

(Tgw)(0  =fxCeAit-T)Bu(r)dr,  t>  0  (126) 

The  Hankel  operator  TG  maps  past  inputs  to  future  outputs,  and  has  finite  rank  n,  that  is,  its 
range  has  finite  dimension  n  if  the  state  space  has  dimension  n  [18],  The  Hankel  singular  values 
are  the  singular  values  of  rG  . 

The  induced  nonn  of  TG  is  defined  by 

|rG|  :=  sup  |rGw|2  (127) 

mgZ-2(-oo,0] 

IMI2-1 

Hankel  model  reduction  is  concerned  with  finding  a  reduced  model  Gr  such  that  its  associated 
Hankel  operator  fG  has  rank,  say  r  <  n ,  in  other  words 


jur:=  min  sup  |Tgw  — rGw| 

r<n  ugL2(-  00, 0]  *  2 

H*1 

=  min  Tg  -Tg 

r<n  '■  r  II 

We  need  the  balanced  realization  for  (73).  Partition  the  state  matrices  as  [54] 

c=[c,  c2] 

21  ^22  L^2  J 

n-rxn-r  _ 

and  the  balanced  empirical  gramians  as 


(128) 


II 

A\  j  Au 

,  B  = 

1 

1 _ 

A  A 

^21  ^22 

-  ^ 
cq 

_  n-rxn-r_ 

(129) 
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rxr 


n-rxn-r . 


The  Hankel  reduced  order  model  has  then  state  space  realization  [54]  [18] 


(130) 


x(t)  =  Arx(t)  +  Bru(t) 
y(t)  =  Crx(t)  +  Dru(t) 

where 

A  =r-V,.+124  +s124A -^+1c;ub;) 
Br=r-\x12Bl+*r+lc;u) 

Cr=C^n+ar+lUB; 

D=-crr+xU 

And  ar+l  is  the  r  + 1  th  Hankel  singular  value  and  U  is  a  unitary  matrix  such  that  [18] 

b2  =  -c\u 
r=z112-^1/ 

The  upper  bound  in  terms  of  the  induced  norm  is  given  by 


(131) 

(132) 

(133) 

(134) 

(135) 

(136) 

(137) 


||G-Gr||:=  sup  \\Gu-Gru\\ 

ugL2  (-oo, oo ) 

HlA1 

<  (o-r+i+o-r+2+  ...  +  o-„) 

In  the  next  section,  this  procedure  is  applied  to  the  nonlinear  convective  flow. 


(138) 


IX,  Application  to  the  Prototype  Nonlinear  Convective  Flow 

The  empirical  balanced  truncation  based  on  linear  systems  is  applied  to  the  Galerkin  model 

a  =  Aa  +  Bu  +  N(a)  +  F,  a(0)  =  ao  (139) 

which  has  an  equilibrium  in  steady  state,  denoted  by  ass . 


The  results  obtained  for  the  first  eight  temporal  coefficients  for  the  full  order  and  the  8  order 
Hankel  reduced  model  are  shown  in  Figure  9.  The  Figure  shows  excellent  agreement  between 
them.  In  Figure  10,  we  compare  the  Hankel  singular  values  of  the  full  order  linearized  and 
Hankel  reduced  8th  order  model.  In  Figure  11,  the  error  between  singular  values  of  linearized 
full  order  and  Hankel  reduced  models  are  plotted. 

Note  that  from  Figure  9  the  reduced  order  Hankel  model  captures  the  largest  singular  values 
which  matter  most  for  the  input-output  behavior  and  control  more  accurately.  The  larger  the 
order  of  the  Hankel  reduced  model  the  smaller  the  error. 
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In  Figure  12  we  compare  the  full  order  solution  of  the  Burgers’  equation  with  the  solution  based 
on  the  8th  order  Flankel  reduced  model.  Figure  12  shows  that  they  behave  similarly  especially  at 
the  boundary  where  control  is  applied. 

In  Figure  13,  we  compare  the  full  order  and  reduced  order  solutions  for  both  the  empirical  and 
Hankel  models.  The  Figure  shows  that  the  lowest  8th  order  Flankel  model  retains  the  important 
behaviors  of  the  flow. 


Figure  9.  Projected  and  POD  Model  Coefficients. 
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Figure  10.  Comparison  of  Hankel  Singular  Values. 
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Figure  11.  Error  Between  Singular  values. 
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Figure  12.  Full  and  Hankel  Reduced  Order  Models  Flow  Responses 
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Figure  13.  Full  and  Reduced  Order  Models’  Responses. 


X.  hT  Controller  Design 

In  this  section  we  consider  the  design  of  an  Hm  controller  using  the  Hankel  reduced  model.  The 
motivation  behind  our  choice  is  that  H°°  controllers  are  robust  against  unmodeled  or  neglected 
dynamics,  and  unknown  or  unmeasurable  disturbances  [18]  [19].  A  fixed  reference  signal 
wref  (x,y)  is  specified  for  the  full  order  system.  Projecting  wref  (x,y)  onto  the  POD  basis  yields 

tracking  coefficients  for  the  reduced  order  model,  denoted  by  aref  .  The  tracking  problem  is 
depicted  in  Figure  16,  where  C  is  the  controller  and  P  the  plant  represented  by  the  dynamical 
equation  (122).  The  computation  of  the  H°°  is  based  on  the  8th  order  Hankel  reduced  model 
(131). 


— o - * 


c 


Figure  14.  Closed-Loop  Control  System  for  Tracking 
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From  Figure  14,  for  tracking  purposes  the  controlled  output  is  chosen  to  be  the  error  signal  e 
which  is  defined  to  be  the  difference  between  the  reference  aref  and  the  actual  output  y(t),  i.e., 

e(t)\=a(t)-aref  (140) 

The  dynamics  of  the  reduced  model  for  tracking  control  represented  in  Figure  16  are  then  given 
by  the  state  space  equation 


x(t)  =  Arx(t )  +  B^aref  +  B2u(t ) 

y(t)  =  Crx(t )  +  Dnaref  +  Duu(t )  (141) 

e(t)  =  - Crx(t )  +  D2lareff  +  D22u{t ) 

The  objective  of  the  II ''  controller  C  is  to  stabilize  the  closed-loop  system  and  minimize  the 
effect  of  aref  on  the  error  e  by  viewing  aref  as  an  unknown  disturbance  in  L2  of  unit  ||  •  ||  - 

norm.  From  Figure  14,  in  terms  of  transfer  function  matrices  of  P  and  C,  the  transfer  matrix 
from  aref  to  e  is  given  by  the  sensitivity  function  Tea  defined  by 

Te^:=(I  +  PCylaref  (142) 

We  compute  the  worst-case  disturbance  transmission  error  due  to  aref  ,  i.e., 

sup  ||e||2  (143) 


which  is  given  by 


sup  e  =ess  sup  a 


0<CO<OO 


T 


(144) 


where  ess  sup  denotes  the  essential  supremum,  and  <r(*)  the  maximum  singular  value  of  its 
argument. 

The  Hx  control  design  reduces  to  the  following  optimization:  Find  C  such  that  the  closed-loop 
system  is  robustly  stable  and 


//:=min 

c 


(145) 


The  solution  of  (145)  is  textbook  material.  There  are  Riccati-based  and  linear  matrix  inequalities 
(LMIs)  based  techniques  to  solve  (1 13)  [18].  In  this  work,  we  use  the  LMI  approach  proposed  in 
Matlab  LMI  toolbox  because  of  its  numerical  robustness  and  stability. 
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After  computing  the  H™  controller  we  close  the  loop  on  the  original  full  order  model.  The 
controller  is  only  8th  order  since  based  on  the  8th  order  Hankel  reduced  model.  The  controlled 
closed-loop  temporal  coefficients  { ak }  are  plotted  for  the  first  8  coefficients  in  Figure  15. 
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Figure  15.  Full  Order  Closed-Loop  System  Tracking  Response 

The  projected  onto  the  POD  basis  initial  condition  and  reference  shown  in  Figure  16. 
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Figure  16.  Flow  Initial  Condition  and  Reference 
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The  controlled  flow  with  the  action  of  the  boundary  controller  is  shown  in  Figure  17,  which 
shows  that  the  flow  tracks  the  reference. 


Figure  17.  Controlled  Flow 


XI.  Model  Estimation  for  NACA  4412  POD  Coefficients  Based  on  Raw  Pressure 
Surface  Sensors  Measurements 


In  this  section,  we  consider  a  2D  flow  control  problem  over  the  NACA  4412  airfoil  [55] [56] [57]. 
The  latter  is  a  model  constructed  from  PIV  experimental  measurements.  The  motivation  for  this 
problem  is  minimization  of  aero-optic  distortion  in  a  bluff-body  flow.  Flere  we  are  primarily 
concerned  with  model  estimation  of  the  temporal  POD  coefficients  using  raw  experimental 
measurements  provided  by  Glauser’s  flow  control  group  at  Syracuse  University.  The  POD 
coefficients  data  are  obtained  from  stochastic  estimation  based  on  1 1  pressure  sensors  along  the 
chord  (surface  measurements)  as  shown  in  Figure  18  [57]. 
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Figure  18.  NACA  4412  Airfoil 


The  NACA  4412  experiment  parameters  are  [57] 
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•  The  Flow  speed:  Uco  =  10  m.s 

•  Airfoil  chord:  20  cm 

•  Reynolds  number:  Re  ~  135,000 

•  Natural  stall  angle:  a  =  15° 


We  describe  the  procedure  employed  to  estimate  the  POD  model  parameters  and  states 
associated  with  the  model  using  the  EM  algorithm  [58]  together  with  Kalman  filtering  [59].  We 
use  the  state  space  model 


Xt+1  =  AtX,+B,Wt 

yt=Ctx,+D,vt 


(146) 


where  xt  e  93"  is  a  state  vector,  yt  e  93"  is  a  measurement  vector,  w,  e  93"'  is  a  state  noise,  and 
v,  €  93rf  is  a  measurement  noise.  The  noise  processes  wt  and  vt  are  assumed  to  be  independent 
zero  mean  and  unit  variance  Gaussian  processes.  Further,  the  noises  are  independent  of  the  initial 
state  x0 ,  which  is  assumed  to  be  Gaussian  distributed. 

The  unknown  system  parameters  Qt  =  {At,Bt,Ct,Dt}  as  well  as  the  system  states  xt  are  estimated 

through  the  measurement  data,  YN  =  { y|,y2,...,yv}  obtained  from  surface  pressure  sensors.  The 

methodology  employed  is  recursive  and  based  on  the  EM  algorithm  together  with  the  Kalman 
filter.  This  is  particularly  useful  for  on-line  implementation  with  the  objective  of  real  time 
control.  The  Kalman  filter  is  discussed  next. 

A.  Model  State  Estimation:  The  Kalman  Filter 

The  Kalman  filter  estimates  the  channel  states  xt  for  given  system  parameter  0  =  {A,B,C,D} 
and  measurements  YN.  It  is  described  by  the  following  equations  [59] 

4  =  + p,\tcTD2  (yt  -  CAxt-x\t-i ) 

xt\t-\  =  Axt-\\t-\  (147) 

A)|o  =  mo 

where  t  =  0, 1, 2,  •  •  • ,  N ,  and  Rt  is  given  by: 

Pti  =  P^+ATB-2A 

P~l  =  CtD  2C  +  B  2  -  B-2Pt]tArB  2  (148) 

Pf\t-\  =  APt_\\,_\AT  +  B2 
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The  channel  parameters  0  =  [A,B,C,D]  are  estimated  using  the  EM  algorithm  which  is 
introduced  next. 

B.  Model  Parameter  Estimation:  The  EM  Algorithm 

The  EM  algorithm  uses  a  bank  of  Kalman  filters  to  yield  a  maximum  likelihood  (ML)  parameter 
estimate  of  the  Gaussian  state  space  model. 

Let  P0  denote  a  fixed  probability  measure;  and  { Pe;  Oe®}  denotes  a  family  of  probability 
measures  induced  by  the  system  parameters  0  .  If  the  original  model  is  a  white  noise  sequence, 
then  { P0;  0  e  0 }  is  absolutely  continuous  with  respect  to  P0  [60].  Moreover,  it  can  be  shown 

that  under  P0  we  have: 


(149) 


The  EM  algorithm  is  an  iterative  scheme  for  computing  the  ML  estimate  of  the  system 
parameters  6 ,  given  the  data  YN .  Specifically,  each  iteration  of  the  EM  algorithm  consists  of  two 
steps:  The  expectation  step  and  the  maximization  step. 

The  expectation  step  evaluates  the  conditional  expectation  of  the  log-likelihood  function  given 
the  complete  data,  and  is  described  by 


A($>  &i)  =  Eg,  | log 1  Yn |  (150) 

where  0,  denotes  the  estimated  system  parameters  at  the  /th  iteration.  The  maximization  step 
finds: 

0l+]  eargmax  A.(G,  0,)  (151) 

The  expectation  and  maximization  steps  are  repeated  until  the  sequence  of  model  parameters 
converge  by  testing  \\0l  -  0/+1||  to  be  less  than  the  required  accuracy. 

The  EM  algorithm  is  described  by  the  following  equations  [58] 
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(152) 


where  B2  =  BBT ,  /J2  =  DDT ,  E(-)  denotes  the  expectation  operator.  These  system  parameters 
A,B2,C,D2 1  can  be  computed  from  the  following  conditional  expectations  [58] 

cAix0*'17"} 

=4z*«0x«iy4 

^ =4z[^fo«+Y,*\]inr} 

C  =  £{zh>,+/s^,]|yj 


where  f?  and  5’  are  given  by: 


Q  = 


T  T 

eiej  +ejei 


;  ij  =  1,2,. ..n 


T 

A,  . 


R  =  s-r-;  i,j  =  1,2,. ..n 


e.e 


S  =  <  -L-s-;  i  =  1,2, ...nr,  n  =  1,2 ,..d 


(153) 


(154) 


in  which  ei  is  the  unit  vector  in  the  Euclidean  space;  that  is  e  =  1  in  the  zth  position,  and  0 
elsewhere.  For  instance,  consider  the  case  n  =  m  =  2,  then  e[ 


(  N 

X 

14 

\t= 1 

J 
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E\ZvL\r„  = 


£!?«,)  £!?(*,*) 

)  C(^) 


(155) 


e,e. 

where  f?,y  =  ■  — ^— ;  i,j  =  1,2| .  The  other  terms  in  (155)  can  be  computed  similarly. 

The  conditional  expectations  j L1'/, L1^’, iff, L'(J’ j  can  be  estimated  from  measurements  YN  as 
follows: 


1)  Filter  estimate  of  ff  : 


L^=E\±x:Qxt\YN 


=  -  )  Tr  (NfPm )  -  i  f  Tr  ( ) 

/  Z  ,=1 

i  (v  r-2rrT,_V<1)  +2xt  P^ra)-xrNa)x 


2^{+xLlB-2ARtN^ATB^2xt]t_1 


(156) 


where  Tr(-)  denotes  the  matrix  trace.  In  (1 15),  r/"  and  N]X)  satisfy  the  following  recursions: 

V"  =  A-P„CTD-2CA)r;»  +  2P„Qx,t_] 

-p,lX'pltcTD'2(y/-cx«,->) 


2)  Filter  estimate  of  L\ ^  : 


=  0„, 


n;u  =  b-\4p,X-:p,,atb-2-2Q 


N(1)  =  0 

Y  0  vmxm 


(157) 


L^=E\YjxTt_lQxt_x\YN 


=  Ee  {40*0  I  Yn}+Eb  j  IX-i Qx,-\  I  Yt 


=  Eg  {xT0  Qx^  +  EA^  x^Qx,  \YN\-Ee  {xTNQxN  \  YN  } 


(158) 


Therefore,  Bl 1  can  be  obtained  from  iff . 
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3)  Filter  estimate  of  L(2)  : 


z®=£{X(^.,+^A)  l^j 

^  ^  t= 1 


1  N 

is 

z  t= 1 


In  this  case,  r/3)  and  satisfy  the  following  recursions: 

V/3)  =  (^  -  -  PAtN?PtVCTD-2  (yt  -  CxtM ) 

+(2P.ltR+2  V24^)h+1 

<43>=a,3> 

/o<3)  =omxl 


'n^=B-2AP,X-IP^TB-2  -2RPtllATB  2 
<  -2B-2APtVRT 

K3)  =  omxm 


(159) 


(160) 


4)  Filter  estimate  of  /) 4  ’ : 


=£(4W4)-4-.^,) 


(161) 


where  r;<4)  satisfy  the  following  recursions: 

V4>  =  fyf  -  P#Cr^G4)r£>  +  2/^Sy, 

■r^X'  (162) 

r(4)  =  0 

/0  umx  1 

Using  the  filters  for  Zfy  (/  =  1,2, 3,4)  and  the  Kalman  filter  described  earlier,  the  system 

parameters  6t  =  [At,Bt,Ct,Dt]  are  estimated  through  the  EM  algorithm  described  in  (155). 

Numerical  results  that  show  the  applicability  of  the  above  algorithm  in  estimating  the  POD 
coefficients  model  parameters  as  well  as  the  model  states  are  showed  below. 
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The  proposed  estimation  algorithm  was  used  in  conjunction  with  1200  samples  of  raw  data  (non- 
liltcred)  which  are  used  to  estimate  recursively  the  model.  Figure  19  and  Figure  20  show  the  two 
first  POD  coefficients.  It  can  be  seen  that  there  is  excellent  agreement  with  the  outputs  of  the 
proposed  model  and  raw  data.  That  is,  once  the  model  is  in  place;  the  raw  measurement  data  can 
be  generated  using  the  proposed  model.  In  Figure  21  we  zoom  in  on  the  2nd  mode. 
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Figure  19.  First  Mode  Coefficient 
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Figure  20.  Second  Mode  Coefficient 
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Figure  21.  Zoom  in  on  Second  Mode  Coefficient 


XII.  Geometric  Interpretation:  N-Width  Approximation 

Geometric  interpretation  of  POD  and  balanced  truncation  in  terms  of  optimizing  the 
Kolmogorov,  Gel’fand  and  linear  n-  widths  of  the  corresponding  compact  operators  is  discussed. 
These  n -widths  quantify  the  inherent  error  generated  in  the  information  collecting  stage  of 
simulation  or  identification,  due  to  lack  of  data  and  inaccurate  measurements,  and  the 
representation  error  due  to  the  loss  of  information  in  the  information  processing  stage. 

The  eigenvalues  ,  ’s  of  (T*Ty2  (or  singular  values  of  T)  defined  in  (13),  and  the  Hankel 
singular  values  <7j  ’s  of  Tc  have  a  geometric  interpretation  in  terms  of  the  computation  of  the  n- 
widths  of  compact  operators  T  and  T0  that  are  defined  on  Hilbert  spaces  Z/([0,  T))  and 
Z2(-oo,0],  respectively.  In  this  section,  we  discuss  the  role  of  POD  and  balanced  truncation  in 
optimizing  different  //-widths  defined  in  [61]  (and  references  therein.) 

In  model  reduction  we  are  interested  in  approximation  by  finite  dimensional  models,  and  in 
particular,  //-parameter  affine  models  such  as  in  POD  and  balanced  truncation.  This  is  related  to 
the  Kolmogorov  /7-width  of  T(L2(Q))  into  Z/([0,  T))  as  the  optimization  which  characterizes  the 
representation  [61] 
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where  X n  is  an  //-dimensional  subspace  of  72([0,  7 )) . 

The  Kolmogorov  n-width  measures  the  extent  to  which  the  space  72([0,7))  can  be  approximated 
by  //-dimensional  subspaces  of  7(7 (O)) ,  it  is  a  measure  of  the  “massivity”  of  7(72(Q)).  It 
represents  the  minimum  representation  error  of  7(72(Q))  by  the  //-dimensional  subspace  Xn  of 
72([0,  7 )) .  In  other  words,  the  Kolmogorov  //-width  quantifies  the  representation  error  due  to 
inaccurate  representation  of  the  set  7(T2(Q)):  It  represents  the  loss  of  information  in  the 
information  processing  stage.  The  inverse  function  of  dn(m)  was  called  the  metric  dimension 
function  by  Zames  in  1976,  and  viewed  as  an  appropriate  measure  of  the  metric  complexity  of 
uncertain  systems.  In  our  case  it  is  the  dimension  of  the  smallest  subspace  whose  elements  can 
approximate  arbitrary  vectors  of  7(72( Q))  to  a  specified  tolerance. 

The  //-width  in  the  sense  of  GeTfand,  is  defined  as 

^(7(Z2(Q));Z2([0,7)))  =  inf  sup  ||  Tf\\2  (164) 

L"  II/Hj  <l,/e£” 

where  the  infimum  is  taken  over  all  subspaces  7"  of  7(Z2(Q))  of  codimension  at  most  //.  If 
J"(7(72(Q));72([0,7))  =  sup{||  /  ||2:  /  e  7(Z2(Q))nZ"} 
where  7'  is  a  subspace  of  codimension  at  most  n,  then  7'  is  an  optimal  subspace  for 

d”  (7  (Z2  (Cl));  Z2  ([0,  7)) 

A  subspace  7"  is  of  codimension  //  if  there  exist  //  continuous  linear  functionals  on 

72  ([0,  7))  for  which 

7"  =  {g  :  g  e  72  (//),  ft  (g)  =  0,  /  =  1, 2, . . . ,  n}  (165) 

The  GeTfand  //-width  characterizes  the  experimental  complexity  of  the  information  collecting 
stage  using  simulation  or  identification.  It  is  related  to  the  inherent  error  due  to  lack  of  data  and 
inaccurate  measurements.  The  inverse  of  the  GeTfand  n-width  gives  the  least  number  of 
measurements  needed  to  reduce  the  modeling  uncertainty  to  a  predetennined  value. 

The  linear  //-width  is  defined  is  defined  by 

^(7(72(Q));72((-^,0])  :=  inf  sup  ||  T0-PJ  ||2 

P«  ||(!>||2<l,(!>el2(0) 

where  P  is  any  continuous  linear  operator  from  f2(Q)  into  f2([0, 7))  of  rank  at  most  //. 

Remark:  Note  similar  definitions  hold  for  the  Hankel  operator  range  TG(72[-oo,  0)) . 

The  basic  results  of  this  section  are  the  following  theorems  which  tell  us  that  the  different  // 
widths  can  be  computed,  and  provide  us  with  explicit  optimal  subspaces  and  operators. 

Theorem:  Let  the  operator  7  be  defined  as  above,  and  let  {/L} ,  {af} ,  {(p.)  be  defined  as  above. 
Then 

dn  (7(72  (Q));  72  ((-<x>,  0]))  =  d "  (7(72  (Q));  72  ((-«,  0]))  =  Sn  (7(72  (Q));  72  ((-«,  0]))  = 

n  =0,  1,  2,--- 
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Furthermore,  the  temporal  coefficients  {«,  }  and  POD  basis  { (pt }  are  optimal  for  the  //-widths  in 
the  following  sense: 

i)  The  subspace  spanned  by  the  coefficients  { ai  J ,  is  optimal  for  X  n  =  Span  { a,,-  ■  ■  ,an  J  , 
is  optimal  for  dH(r(Z2(Q));Z2((-oo,0]))  . 

ii)  The  subspace  Ln  =  {</>  e  L2  (Q),  <  (/),  cpt  >j  =  0,  i  =  1, 2,  •  •  • ,  n }  is  optimal  for 

dn(T(L2(Cl));L2((-K,  0])). 

iii)  The  linear  operator  Pn(f)  -  y,.' ,  <  (j),  (p;  >,  (pi  is  optimal  for  £h(T(Z2(Q));Z2((-oo,  0])). 
A  similar  Theorem  holds  for  the  Hankel  operator  TG  and  is  stated  next 

Theorem:  Let  the  operator  TG  be  defined  as  above,  and  let  {cr  }  ,  {x,}  ,{C,}  be  defined  as  above. 
Then 

d„  (rG  (L2  ((-00, 0]));  f2((0,oo]))  =  d n  (rG  {Is  ((-«,  0]));  Is  ((0,  oo])) 

=  8n  (r  G  (L2  ((-oo,  0]));  Z2  ((0,  co])) 

=  T„+1,  «  =0, 1,  2,  — 

F urthermo re ,  the  temporal  coefficients  lx,}  and  POD  basis  { C,  J  are  optimal  for  the  //-widths  in 
the  following  sense: 

i)  The  subspace  spanned  by  the  coefficients  { <£}  ,  Xn  =  Span  { f  },  is  optimal  for 

for  dn  (T G  (i}  ((-oo,  0]));  L2  ((0,  oo])) 

T  f0 

ii)  The  subspace  L'  =  (x  eL"(-oo,0],  X(r)X, (r)  - =  1,2, ■•■,//}  is  optimal  for 

J-oo 

J''(rG(L2((-»,0]));Z2((0,»])). 

i  f° 

iii)  The  linear  operator  Qn(!)  =  /_i  J  </>(T)%i(.T)dTXi  is  optimal  for 

1=1  _C° 

d„  (r G  (X2  (-°°,  0]);  Is  [0,  oo))  . 

POD  may  fail  to  capture  the  nonlinear  degrees  of  freedom  in  nonlinear  PDEs,  in  particular  for 
navier-Stokes  equations  with  Reynolds  numbers  in  the  thousands,  since  it  assumes  that  data 
belong  to  a  linear  space  and  therefore  relies  on  the  Euclidean  distance  as  the  metric  to  minimize. 
However,  snapshots  generated  by  nonlinear  PDEs  belong  to  manifolds  for  which  the  geodesics, 
when  they  exist,  do  not  correspond  in  general  to  the  Euclidean  distance.  A  geodesic  is  a  curve 
that  is  locally  the  shortest  path  between  two  points.  In  this  next  section,  we  propose  a  model 
reduction  method  which  generalizes  POD  to  nonlinear  fluid  flows  corresponding  to  manifolds 
which  have  a  differentiable  structure  at  each  of  their  points.  The  algorithm  is  applied  to  a 
prototype  nonlinear  turbulent  flow  governed  by  the  Navier-Stokes  equation. 

XIII.  Nonlinear  Proper  Orthogonal  Decomposition  (POD):  Auto-Associative  Models 
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Nonlinear  POD  is  inspired  from  auto-associative  models  that  have  been  introduced  as  a  new  tool 
that  deals  with  nonlinear  data.  Such  models  rely  on  successive  approximations  of  a  dataset  by 
manifolds  of  increasing  dimensions.  In  this  paper,  auto-associative  models  are  proposed  as 
candidates  to  generalize  POD  to  nonlinear  systems.  These  models  are  dedicated  to  the 
approximation  of  datasets  by  manifolds  [62],  Let  us  introduce  definitions 


A  function  Fd  :  IRA  — »  RA  is  a  (/-dimensional  auto-associative  function  if  there  exist  d  unit 
orthogonal  vectors  ak ,  called  proper  directions,  and  d  continuously  differentiable  functions 
sk  :  R  — »  !A ,  called  regression  functions,  such  that  [62] 

Pgj  °sk  =  Sj  kIdM  for  all  1  <  j  <k  <d  (166) 

where  5j  k  is  the  Kronecker  symbol  and  for  each  unit  vector  a  et'',  we  denote  by  P  (•)  =  (a,  •) 


the  linear  projection  from  1A  to  R  .  Also,  for  all  sets  E  ,  the  identity  function  E  — >  E  is  denoted 
by  Id E  [62].  Define 


Fd  ~  (IdRp  -sd  ° Pa,,)o---° (IdRP  -s1  °Pgl) 

=  11  Wv-sk°Pj) 


(167) 


The  equation  Fd(x)  =  0,  x  defines  a  differentiable  d-dimensional  manifold  of  R''  [63], 

Thus,  the  equation  Fd(x)  =  0  defines  a  space  in  which  every  point  has  a  neighborhood  which 

resembles  the  Euclidean  space  W1 ,  but  in  which  the  global  structure  may  be  more  complicated. 
As  an  example,  on  a  1 -dimensional  manifold,  every  point  has  a  neighborhood  that  resembles  a 
line.  In  a  2-manifold,  every  point  has  a  neighborhood  that  looks  like  a  plane. 

Let  X  be  a  square  integrable  vector  function  of  1L  X  represents  snapshots.  For  an  auto- 
associative  function  Fd ,  considers  =  F"  (X)  where  in  this  context,  ||^||]  is  called  the  residual 

error.  PDE  Snapshots  X  are  approximated  by  the  manifold  Fd(x)  =  0,  xgK?  and  ||^||] 
represents  the  (squared)  norm  of  X  outside  the  manifold. 

Note  that  such  vector  X  always  satisfies  a  0-dimensional  auto-associative  model  with  F()  =  Id  p 
and  l^l]  =  ||x||] .  Similarly,  X  always  satisfies  a  /^-dimensional  auto-associative  model  with 

II  ||2 

Fp  =  0  and  ||£'||  ,  =  0 .  In  practice,  it  is  important  to  find  a  balance  between  these  two  extreme 

II  ||2  II  m2 

cases  by  constructing  a  ^-dimensional  model  with  d  <k  p  and  ||f  |  <sc  \\X\\2 . 

For  instance,  in  the  case  where  the  autocorrelation  X  of  X  is  of  rank  d,  then  X  is  located  on  a 
J-dimensional  linear  subspace  defined  by  the  equation  FdOD(x )  =  0  with 

FPdOD(x)  =  x-fjPk(x)ak  (168) 

k= 1 
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where  cik ,k  =  \,---,d  are  the  eigenvectors  of  X  associated  to  the  positive  eigenvalues.  Equation 
(168)  can  be  rewritten  as  Fd (x)  =  0,  where  Fd is  a  (/-dimensional  auto-associative  function  with 
linear  regression  functions  ak  ( t ),  k  =  1, — ,  d .  Moreover,  we  have  ||£j|]  =  0 . 

Since  (168)  is  the  model  produced  by  a  POD,  it  straightforwardly  follows  that  POD  is  a  special 
(linear)  case  of  auto-associative  models.  In  the  next  section,  we  show  an  algorithm  to  build  auto- 
associative  models  with  non-necessarily  linear  regression  functions,  small  dimension  and  small 
residual  error. 

A.  Projection  Algorithm 

Given  a  unit  vector  aeMp,  an  index  7:M— »Mis  a  functional  measure  the  interest  of  the 
projection  Pa(X)  with  a  non  negative  real  number.  This  depends  on  the  considered  data  analysis 

problem.  For  instance,  a  possible  choice  of  /  is  the  projected  squared  L2  nonn  I°Pa(»)  =  ||Pa(»)||2 . 
So,  the  maximization  of  I °Pa(X)  with  respect  to  a  yields  the  most  interesting  direction  for  this 
given  criteria  [64]. 

Let  d  e  {0,  ••■,/?},  the  projection  algorithm  applies  the  following  iteratively:  computation  of  the 
axes,  projection,  regression  and  update  [64][65]. 

Algorithm:  Define  R°  =  X 
For  k  =  !,■■■,  d  : 

1)  Determine  ak  -  argmax/  °P  (Rk~l)  such  that  llxll  =  l,P2  (x)  =  0,1  <  j  <  k 

2)  Compute  Yk  =  Pt  (Rk  1 ) 

3)  Estimate  sk  (t)Yk  /  =  Rk  1 

4)  Compute  Rk  =Rk~l  -sk(Yk) 

Step  1  computes  an  axis  orthogonal  to  the  previous  ones  and  maximizes  a  given  index/  .  Step  2 
projects  the  residuals  on  this  axis  to  determine  the  proper  variables  Yk ,  and  step  3  estimates  of 
the  regression  function  of  FA  best  approximating  the  residuals.  Finally  Step  4  updates  the 
residuals. 

In  the  next  section,  the  algorithm  is  applied  to  a  prototype  nonlinear  turbulent  flow  governed  by 
the  Navier-Stokes  equation. 

B.  Full  Order  System 

To  test  the  algorithm,  our  data  is  a  numerical  finite  element  solution  of  Navier  Stokes  equations 
that  describes  the  dynamics  of  an  incompressible  turbulent  flow.  General  Navier  -Stokes 
equations  are  of  the  following  form  [66] 
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(169) 


dtp  +  V-\pv\  =  0, 


juVv  +—juV -vl 


d,(Pv)  +  pv-Vv-V 
d,  ( cpPT )  +  cppvVT  -  V  •  [AVT\  =  h 


1 


+  VPto,  =pf 


(170) 

(171) 


where  v  is  the  velocity,  p  the  density,  ptot  the  total  pressure,  and  T  the  temperature  of  the  fluid 
occupying  a  two-  or  three-dimensional  region  .  The  parameters  p  >  0  (dynamic  viscosity), 
cp  >  0  (heat  capacity)  and  A  >  0  (heat  conductivity)  characterize  the  properties  of  the  fluid.  The 
volume  force  /  and  the  heat  source  h  are  explicitly  given. 


In  temperature-driven  flows,  h  may  implicitly  depend  on  the  temperature  and  further  quantities 
describing  heat  release,  as  for  example  by  chemical  reactions.  For  simplicity,  we  assume  an 
incompressible  flow  (i.e.  p  is  constant),  then  we  have: 

V»v  =  0 

In  this  model,  we  consider  as  the  primal  unknowns  the  velocity  v  ,  the  pressure  p  =  ptot ,  and  the 

temperature  T  .  For  most  parts  of  the  discussion,  the  flow  is  assumed  to  be  isothennal,  so  that 
the  energy  equation  decouples  from  the  momentum  and  continuity  equations,  and  the 
temperature  only  enters  through  the  viscosity  parameter.  The  system  is  closed  by  imposing 
appropriate  initial  and  boundary  conditions  for  the  flow  variables 


vu=v 


=  0, 


rigid 


Vlr,„  =V 


(pdnv  +  pn)\r  =0 


and  corresponding  ones  for  the  temperature,  where  Trigid,  r/Hand  Tout  are  the  rigid  part,  the 
inflow  part  and  the  outflow  part  of  the  boundary  o’Q  ,  respectively. 

Assuming  the  isothermal  case,  i.e.  >c>0  =  1 ,  the  Navier  Stokes  system  can  be  written  as  [66]: 

dtv+vVv-kVv+Vp  = /,  V»v  =  0  (171) 

with  the  kinematic  viscosity  parameter  k  =  —  .  In  this  formulation  the  domain  Q  may  be  taken 

Po 

two  or  three  dimensional  according  to  the  particular  requirements  of  the  simulation. 

C.  Weak  Formulation 

Taking  the  inner  product  of  both  sides  of  (171)  with  the  i-th  nonlinear  POD  mode  Yl  and 
utilizing  Green's  identities  results  in  the  weak  formulation  in  a  similar  fashion  as  in  [12]  to  get: 
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(172) 


(<3,v,  (p)  +  (' Vv,  Wcp)  +  (v.Vv,  <p)~(p,  V.v)  =  (/,  p) 
re  (•,  •)  denote  the  inner  product  in  L2(Q) ,  and  smooth  cp  e  Lr{D.) . 

The  problem  is  numerically  solved  using  Comsol  Multiphysics  [66].  Figure  24  shows  the 
problem  geometry,  boundary  conditions  and  the  full  order  velocity  solution.  Wortmann  FX  60- 
100  airfoil  is  used.  The  finite  element  discretization  results  in  10029  degree  of  freedom,  other 
parameters  are  as  follows: 

Density:  p  =  lkg/m 3 

Dynamic  Viscosity:  p  =  O.OOlPascal.sec 
Reynolds  Number:  Re=1000. 

The  full  order  fluid  flow  model  is  represented  in  Figure  22.  The  reduced  order  model  has  50 
states  and  is  represented  in  Figure  23.  Both  Figures  show  that  good  agreement  between  both 
flows.  We  are  currently  generalizing  the  method  to  Reynolds  numbers  in  the  tens  of  thousands. 


Figure  22.  Full  Order  Finite  Element  Solution 


Figure  23.  Reduced  Order  Model 


XIV.  Conclusion 

In  this  project,  model  reduction  algorithms  and  control  design  approaches  for  aerodynamic  flows 
were  developed.  The  goal  of  this  research  is  eventually  to  develop  a  methodology  for  feedback 
control  of  fluid  flows  modeled  by  Navier-Stokes  equations.  The  proposed  work  places  an 
emphasis  upon  aerodynamic  control.  For  the  case  of  boundary  control,  the  development  of 
reduced  models  has  been  an  open  problem.  For  many  Air  Force  applications  of  practical  interest, 
such  as  feedback  control  of  the  air  flow  over  an  airplane  wing,  boundary  actuation  is  a 
requirement.  These  applications  require  actuation  to  be  located  on  the  surface  if  they  are  to  be 
implemented  in  hardware  in  the  physical  system.  Most  reduced  modeling  efforts  have  either  been 
concerned  with  the  case  of  control  action  via  a  body  force  or  have  approximated  boundary 
actuation  by  a  body  force  on  the  domain  interior.  In  this  project,  we  have  developed  a 
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methodology  which  extracts  boundary  conditions  in  reduced  order  POD  and  finite  difference 
models.  For  the  most  part,  this  method  utilizes  POD  and  finite  difference  bases  in  combination 
with  a  weak  formulation  of  the  Galerkin  projection  to  construct  models  where  boundary  control 
input  appears  explicitly.  Special  care  is  taken  when  POD  is  used  for  the  development  of 
feedback  control  laws.  In  particular,  construction  of  a  POD  basis  capable  of  spanning  the 
baseline  solution,  as  well  as  dynamics  introduced  by  boundary  actuation  has  been  addressed 
using  a  method  called  split  POD. 

Tools  borrowed  from  the  theory  of  operators  were  used  to  show  that  POD  and  balanced 
truncation  are  optimal  in  a  precise  sense.  Optimality  is  quantified  in  terms  of  optimal 
approximations  by  lower  rank  Hilbert-Schmidt  (integral)  operators  in  Hilbert-Schmidt  norms. 
The  difference  in  the  two  model  reduction  techniques  lies  in  the  fact,  that  the  optimizations  occur 
in  different  integral  operators  defined  on  different  Is  spaces.  However,  both  optimal 
approximations  are  posed  in  Hilbert  operator  spaces  allowing  for  the  optimal  approximations  to 
be  computed  explicitly. 

Model  estimation  and  identification  of  the  POD  modes  for  a  turbulent  flow  over  the  NACA  4412 
airfoil  is  developed.  The  model  is  constructed  from  PIV  experimental  measurements  obtained 
from  the  flow  control  group  at  Syracuse  University.  The  dynamics  of  the  flow  are  captured  using 
a  linear  stochastic  state-space  model.  Estimation  and  identification  algorithms  solely  from  wing 
surface  pressure  measurement  data  are  developed.  These  algorithms  are  based  on  combining  the 
KF  with  the  EM  algorithm  that  estimate  and  identify  recursively  the  state  of  the  flow  and  its 
parameters,  respectively.  Numerical  results  are  provided  to  evaluate  the  accuracy  of  the  proposed 
algorithms. 

Empirical  balanced  truncation  and  Hankel  nonn  model  reduction  have  been  considered  in 
conjunction  with  POD  as  an  approach  for  deriving  reduced-order  models.  Like  POD,  both  are 
based  on  simulation/experimental  data  and  can  be  implemented  via  standard  matrix 
computations.  Improvements  to  the  scheme  originally  proposed  in  [5]  [23]  have  been  presented 
that  lead  to  reduced  data  requirements  that  may  become  significant  for  applications  such  as 
aerodynamic  flow  control.  Essentially,  the  balancing  transformation  is  constructed  via  Markov 
parameters  that  can  be  identified  from  measurements  collected  in  a  single  experiment/ 
simulation.  The  approach  has  been  applied  with  favorable  results  to  nonlinear  convection 
governed  by  the  2D  Burgers’  equation,  a  partial  differential  equation  in  two  spatial  dimensions 
that  possesses  features  comparable  to  the  Navier-Stokes  equations  governing  fluid  flow.  A  HJ 
feedback  flow  controller  was  designed  based  on  the  empirical  reduced  model  to  achieve  flow 
tracking.  The  closed-loop  on  the  full  order  model  shows  good  flow  tracking  performance. 

Geometric  interpretation  of  POD  and  balanced  truncation  in  terms  of  optimizing  the 
Kolmogorov,  GeTfand  and  linear  n- widths  is  discussed.  These  «- widths  quantify  the  inherent 
and  representation  errors  generated  in  the  infonnation  collecting  and  processing  stages  in 
simulation  or  identification.  To  the  best  of  our  knowledge  this  is  the  first  time  that  model 
reduction  and  metric  complexity  theory  with  the  notion  of  n-width  are  related  to  each  other  in  an 
explicit  way. 
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A  novel  model  reduction  method  which  generalizes  POD  to  nonlinear  PDEs  with  solutions 
belonging  to  manifolds  which  have  (approximately)  a  differentiable  structure  at  each  of  their 
points.  Conventional  POD  is  widely  used  to  reduce  the  order  of  high  order  models  of  nonlinear 
system  but  it  fails  with  systems  of  high  nonlinearity  because  the  projection  is  done  to  a  linear 
Euclidean  space.  Non-linear  POD  proposed  here  solves  this  limitation  using  a  projection  on  the 
manifold  constructed  from  the  input  data.  In  future  work  we  plan  to  improve  and  test  our 
algorithm  on  systems  with  higher  nonlinearities  by  increasing  the  Reynolds  numbers  of  turbulent 
flows. 
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